Pricing Bermudan options under local Lévy models with default

Anastasia Borovykh Dipartimento di Matematica, Università di Bologna, Bologna, Italy. e-mail: anastasia.borovykh2@unibo.it Andrea Pascucci Dipartimento di Matematica, Università di Bologna, Bologna, Italy. e-mail: andrea.pascucci@unibo.it Cornelis W. Oosterlee Centrum Wiskunde & Informatica, Amsterdam, The Netherlands. e-mail: c.w.oosterlee@cwi.nlDelft University of Technology, Delft, The Netherlands.

(This version: August 7, 2025)

Abstract

We consider a defaultable asset whose risk-neutral pricing dynamics are described by an exponential Lévy-type martingale. This class of models allows for a local volatility, local default intensity and a locally dependent Lévy measure. We present a pricing method for Bermudan options based on an analytical approximation of the characteristic function combined with the COS method. Due to a special form of the obtained characteristic function the price can be computed using a Fast Fourier Transform-based algorithm resulting in a fast and accurate calculation. The Greeks can be computed at almost no additional computational cost. Error bounds for the approximation of the characteristic function as well as for the total option price are given.

Keywords: Bermudan option, local Lévy model, defaultable asset, asymptotic expansion, Fourier-cosine expansion

1 Introduction

In financial mathematics, the fast and accurate pricing of financial derivatives is an important branch of research. Depending on the type of financial derivative, the mathematical task is essentially the computation of integrals, and this sometimes needs to be performed in a recursive way in a time-wise direction. For many stochastic processes that model the financial assets, these integrals can be most efficiently computed in the Fourier domain. However, for some relevant and recent stochastic models the Fourier domain computations are not at all straightforward, as these computations rely on the availability of the characteristic function of the stochastic process (read: the Fourier transform of the transitional probability distribution), which is not known. This is especially true for state-dependent asset price processes, and for asset processes that include the notion of default in their definition. With the derivations and techniques in the present paper we make available the highly efficient pricing of so-called Bermudan options to the above mentioned classes of state-dependent asset dynamics, including jumps in asset prices and the possibility of default. In this sense, the class of asset models for which Fourier option pricing is highly efficient increases by the contents of the present paper. Essentially, we approximate the characteristic function by an advanced Taylor-based expansion in such a way that the resulting characteristic function exhibits favorable properties for the pricing methods.

Fourier methods have often been among the winners in option pricing competitions such as BENCHOP [16]. In [5], a Fourier method called the COS method, as introduced in [4], was extended to the pricing of Bermudan options. The computational efficiency of the method was based on a specific structure of the characteristic function allowing to use the fast Fourier transform (FFT) for calculating the continuation value of the option. Fourier methods can readily be applied to solving problems under asset price dynamics for which the characteristic function is available. This is the case for exponential Lévy models, such as the Merton model developed in [13], the Variance-Gamma model developed in [12], but also for the Heston model [6]. However, in the case of local volatility, default and state-dependent jump measures there is no closed form characteristic function available and the COS method can not be readily applied.

Recently, in [14] the so-called adjoint expansion method for the approximation of the characteristic function in local Lévy models is presented. This method is worked out in the Fourier space by considering the adjoint formulation of the pricing problem, that is using a backward parametrix expansion as was also later done in [1]. In this paper we generalize this method to include a defaultable asset whose risk-neutral pricing dynamics are described by an exponential Lévy-type martingale with a state-dependent jump measure, as has also been considered in [11] and in [7].

Having obtained the analytical approximation for the characteristic function we combine this with the COS method for Bermudan options. We show that this analytical formula for the characteristic function still possesses a structure that allows the use of a FFT-based method in order to calculate the continuation value. This results in an efficient and accurate computation of the Bermudan option value and of the Greeks. The characteristic function approximation used in the COS method is already very accurate for the $2$ nd-order approximation, meaning that the explicit formulas are simple and this makes method easy and quick to implement. Finally, we present a theoretical justification of the accurate performance of the method by giving the error bounds for the approximated characteristic function.

The rest of this paper is organized as follows. In Section 2 we present the general framework which includes a local default intensity, a state-dependent jump measure and a local volatility function. Then we derive the adjoint expansion of the characteristic function. In Section 3 we propose an efficient algorithm for calculating the Bermudan option value, which makes use of the Fast Fourier transform. In Section 4 we prove error bounds for the $0$ th- and $1$ st-order approximation, justifying the accuracy of the method. Finally, in Section 5 numerical examples are presented, showing the flexibility, accuracy and speed of the method.

2 General framework

We consider a defaultable asset $S$ whose risk-neutral dynamics are given by:

$\displaystyle S_{t}$	$\displaystyle={\mathds{1}}_{\{t<\zeta\}}e^{X_{t}},$
$\displaystyle dX_{t}$	$\displaystyle=\mu(t,X_{t})dt+\sigma(t,X_{t})dW_{t}+\int_{\mathbb{R}}d\tilde{N}_{t}(t,X_{t-},dz)z,$
$\displaystyle d\tilde{N}_{t}(t,X_{t-},dz)$	$\displaystyle=dN_{t}(t,X_{t-},dz)-\nu(t,X_{t-},dz)dt,$
$\displaystyle\zeta$	$\displaystyle=\inf\{t\geq 0:\int_{0}^{t}\gamma(s,X_{s})ds\geq\varepsilon\},$	(2.1)

where $\tilde{N}_{t}(t,x,dz)$ is a compensated random measure with state-dependent Lévy measure $\nu(t,x,dz)$ . The default time $\zeta$ of $S$ is defined in a canonical way as the first arrival time of a doubly stochastic Poisson process with local intensity function $\gamma(t,x)\geq 0$ , and $\varepsilon\sim\mathrm{Exp}(1)$ and is independent of $X$ . Thus the model features:

•

a local volatility function $\sigma(t,x)$ ;
•

a local Lévy measure: jumps in $X$ arrive with a state-dependent intensity described by the local Lévy measure $\nu(t,x,dz)$ . The jump intensity and jump distribution can thus change depending on the value of $x$ . A state-dependent Lévy measure is an important feature because it allows to incorporate stochastic jump-intensity into the modeling framework;
•

a local default intensity ${\gamma}(t,x)$ : the asset $S$ can default with a state-dependent default intensity.

This way of modeling default is also considered in a diffusive setting in [3] and for exponential Lévy models in [2].

We define the filtration of the market observer to be $\mathcal{G}=\mathcal{F}^{X}\vee\mathcal{F}^{D}$ , where $\mathcal{F}^{X}$ is the filtration generated by $X$ and $\mathcal{F}_{t}^{D}:=\sigma(\{\zeta\leq u\},u\leq t)$ , for $t\geq 0$ , is the filtration of the default. We assume

\int_{\mathbb{R}}e^{|z|}\nu(t,x,dz)<\infty,

(2.2)

and by imposing that the discounted asset price $\tilde{S}_{t}:=e^{-rt}S_{t}$ is a $\mathcal{G}$ -martingale, we get the following restriction on the drift coefficient:

\mu(t,x)=\gamma(t,x)+r-\frac{\sigma^{2}(t,x)}{2}-\int_{\mathbb{R}}\nu(t,x,dz)(e^{z}-1-z).

Is it well-known (see, for instance, [8, Section 2.2]) that the price $V$ of a European option with maturity $T$ and payoff $\Phi(S_{T})$ is given by

\displaystyle V_{t}={\mathds{1}}_{\{\zeta>t\}}e^{-r(T-t)}E\left[e^{-\int_{t}^{T}\gamma(s,X_{s})ds}{\varphi}(X_{T})|X_{t}\right],\qquad t\leq T,

(2.3)

where ${\varphi}(x)=\Phi(e^{x})$ . Thus, in order to compute the price of an option, we must evaluate functions of the form

\displaystyle u(t,x):=E\left[e^{-\int_{t}^{T}\gamma(s,X_{s})ds}{\varphi}(X_{T})|X_{t}=x\right].

(2.4)

Under standard assumptions, $u$ can be expressed as the classical solution of the following Cauchy problem

\displaystyle\begin{cases}Lu(t,x)=0,\qquad&t\in[0,T[,\ x\in\mathbb{R},\\ u(T,x)={\varphi}(x),&x\in\mathbb{R},\end{cases}

(2.5)

where $L$ is the integro-differential operator

	$\displaystyle Lu(t,x)=$	$\displaystyle\ \partial_{t}u(t,x)+r\partial_{x}u(t,x)+\gamma(t,x)(\partial_{x}u(t,x)-u(t,x))+\frac{\sigma^{2}(t,x)}{2}(\partial_{xx}-\partial_{x})u(t,x)$
		$\displaystyle-\int_{\mathbb{R}}\nu(t,x,dz)(e^{z}-1-z)\partial_{x}u(t,x)+\int_{\mathbb{R}}\nu(t,x,dz)(u(t,x+z)-u(t,x)-z\partial_{x}u(t,x)).$		(2.6)

The function $u$ in (2.4) can be represented as an integral with respect to the transition distribution of the defaultable log-price process $\log S$ :

\displaystyle u(t,x)=\int_{\mathbb{R}}{\varphi}(y)\Gamma(t,x;T,dy).

(2.7)

Here we notice explicitly that $\Gamma(t,x;T,dy)$ is not necessarily a standard probability measure because its integral over ${\mathbb{R}}$ can be strictly less than one; nevertheless, with a slight abuse of notation, we say that its Fourier transform

\hat{\Gamma}(t,x;T,{\xi}):=\mathcal{F}(\Gamma(t,x;T,\cdot))(\xi):=\int_{\mathbb{R}}e^{i\xi y}\Gamma(t,x;T,dy),\qquad{\xi}\in{\mathbb{R}},

is the characteristic function of $\log S$ .

2.1 Adjoint expansion of the characteristic function

In this section we generalize the results in [14] to our framework and develop an expansion of the coefficients

a(t,x):=\frac{\sigma^{2}(t,x)}{2},\qquad\gamma(t,x),\qquad\nu(t,x,dz),

around some point $\bar{x}$ . The coefficients $a(t,x)$ , $\gamma(t,x)$ and $\nu(t,x,dz)$ are assumed to be continuously differentiable with respect to $x$ up to order $N\in\mathbb{N}$ .

From now on for simplicity we assume that the coefficients are independent of $t$ (see Remark 2.2 for the general case). First we introduce the $n$ th-order approximation of $L$ in (2):

	$\displaystyle L_{n}=$	$\displaystyle\ L_{0}+\sum_{k=1}^{n}\Big{(}(x-\bar{x})^{k}a_{k}(\partial_{xx}-\partial_{x})+(x-\bar{x})^{k}\gamma_{k}\partial_{x}-(x-\bar{x})^{k}\gamma_{k}$
		$\displaystyle-\int_{\mathbb{R}}(x-\bar{x})^{k}\nu_{k}(dz)(e^{z}-1-z)\partial_{x}+\int_{\mathbb{R}}(x-\bar{x})^{k}\nu_{k}(dz)(e^{z\partial_{x}}-1-z\partial_{x})\Big{)},$		(2.8)

where

\displaystyle L_{0}

\displaystyle=\partial_{t}+r\partial_{x}+a_{0}(\partial_{xx}-\partial_{x})+\gamma_{0}\partial_{x}-\gamma_{0}-\int_{\mathbb{R}}\nu_{0}(dz)(e^{z}-1-z)\partial_{x}+\int_{\mathbb{R}}\nu_{0}(dz)(e^{z\partial_{x}}-1-z\partial_{x}),

(2.9)

and

\displaystyle a_{k}=\frac{\partial_{x}^{k}a(\bar{x})}{k!},\qquad\gamma_{k}=\frac{\partial_{x}^{k}\gamma(\bar{x})}{k!},\qquad\nu_{k}(dz)=\frac{\partial_{x}^{k}\nu(\bar{x},dz)}{k!},\qquad\ k\geq 0.

(2.10)

The basepoint $\bar{x}$ is a constant parameter which can be chosen freely. In general the simplest choice is $\bar{x}=x$ (the value of the underlying at initial time $t$ ): we will see that in this case the formulas for the Bermudan option valuation are simplified.

Let us assume for a moment that $L_{0}$ has a fundamental solution $G^{0}(t,x;T,y)$ that is defined as the solution of the Cauchy problem

\begin{cases}L_{0}G^{0}(t,x;T,y)=0\qquad&t\in[0,T[,\ x\in\mathbb{R},\\ G^{0}(T,\cdot;T,y)={\delta}_{y}.\end{cases}

In this case we define the $n$ th-order approximation of $\Gamma$ as

\Gamma^{(n)}(t,x;T,y)=\sum_{k=0}^{n}G^{k}(t,x;T,y),

where, for any $k\geq 1$ and $(T,y)$ , $G^{k}(\cdot,\cdot;T,y)$ is defined recursively through the following Cauchy problem

\begin{cases}L_{0}G^{k}(t,x;T,y)=-\sum\limits_{h=1}^{k}(L_{h}-L_{h-1})G^{k-h}(t,x;T,y)\qquad&t\in[0,T[,\ x\in\mathbb{R},\\ G^{k}(T,x;T,y)=0,&x\in\mathbb{R}.\end{cases}

Notice that

	$\displaystyle L_{h}-L_{h-1}=$	$\displaystyle\ (x-\bar{x})^{h}a_{h}(\partial_{xx}-\partial_{x})+(x-\bar{x})^{h}\gamma_{h}\partial_{x}-(x-\bar{x})^{h}\gamma_{h}$		(2.11)
		$\displaystyle-\int_{\mathbb{R}}(x-\bar{x})^{h}\nu_{h}(dz)(e^{z}-1-z)\partial_{x}+\int_{\mathbb{R}}(x-\bar{x})^{h}\nu_{h}(dz)(e^{z\partial_{x}}-1-z\partial_{x}).$		(2.12)

Correspondingly, the $n$ th-order approximation of the characteristic function $\hat{\Gamma}$ is defined to be

\hat{\Gamma}^{(n)}(t,x;T,\xi)=\sum_{k=0}^{n}\mathcal{F}\left(G^{k}(t,x;T,\cdot)\right)(\xi):=\sum_{k=0}^{n}\hat{G}^{k}(t,x;T,\xi),\qquad{\xi}\in{\mathbb{R}}.

(2.13)

Now we remark that the operator $L$ acts on $(t,x)$ while the characteristic function is a Fourier transform taken with respect to $y$ : in order to take advantage of such a transformation, in the following theorem we characterize $\hat{\Gamma}^{(n)}$ in terms of the Fourier transform of the adjoint operator $\tilde{L}=\tilde{L}^{(T,y)}$ of $L$ , acting on $(T,y)$ .

Theorem 2.1 (Dual formulation).

For any $(t,x)\in]0,T]\times\mathbb{R}$ , the function $G^{0}(t,x;\cdot,\cdot)$ is defined through the following dual Cauchy problem

\begin{cases}\tilde{L}_{0}^{(T,y)}G^{0}(t,x;T,y)=0\qquad&T>t,\ y\in\mathbb{R},\\ G^{0}(T,x;T,\cdot)={\delta}_{x}.\end{cases}

(2.14)

where

\displaystyle\tilde{L}_{0}^{(T,y)}

\displaystyle=-\partial_{T}-r\partial_{y}+a_{0}(\partial_{yy}+\partial_{y})-\gamma_{0}\partial_{y}-\gamma_{0}+\int_{\mathbb{R}}\nu_{0}(dz)(e^{z}-1-z)\partial_{y}+\int_{\mathbb{R}}\bar{\nu}_{0}(dz)(e^{z\partial_{y}}-1-z\partial_{y}).

(2.15)

Moreover, for any $k\geq 1$ , the function $G^{k}(t,x;\cdot,\cdot)$ is defined through the dual Cauchy problem as follows:

\begin{cases}\tilde{L}_{0}^{(T,y)}G^{k}(t,x;T,y)=-\sum\limits_{h=1}^{k}\left(\tilde{L}_{h}^{(T,y)}-\tilde{L}_{h-1}^{(T,y)}\right)G^{k-h}(t,x;T,y)\qquad&T>t,\ y\in\mathbb{R},\\ G^{k}(T,x;T,y)=0&y\in\mathbb{R},\end{cases}

(2.16)

with

$\displaystyle\tilde{L}_{h}^{(T,y)}-\tilde{L}_{h-1}^{(T,y)}=$	$\displaystyle\ a_{h}h(h-1)(y-\bar{x})^{h-2}+a_{h}(y-\bar{x})^{h-1}\left(2h\partial_{y}+(y-\bar{x})(\partial_{yy}+\partial_{y})+h\right)$	(2.17)
	$\displaystyle-\gamma_{h}h(y-\bar{x})^{h-1}-\gamma_{h}(y-\bar{x})^{h}\left(\partial_{y}+1\right)$	(2.18)
	$\displaystyle+\int_{\mathbb{R}}\nu_{h}(dz)(e^{z}-1-z)\left(h(y-\bar{x})^{h-1}+(y-\bar{x})^{h}\partial_{y}\right)$
	$\displaystyle+\int_{\mathbb{R}}\bar{\nu}_{h}(dz)\left((y+z-\bar{x})^{h}e^{z\partial_{y}}-(y-\bar{x})^{h}-z\left(h(y-\bar{x})^{h-1}-(y-\bar{x})^{h}\partial_{y}\right)\right),$	(2.19)

where in defining the adjoint of the operator we use the notation

\displaystyle e^{z\partial_{y}}f(y):=\sum_{n=0}^{\infty}\frac{z^{n}}{n!}\partial_{y}^{n}f(y)=f(y+z).

(2.20)

Notice that the adjoint Cauchy problems (2.14) and (2.16) admit a solution in the Fourier space and can be solved explicitly; in fact, we have

\mathcal{F}\left(\tilde{L}_{0}^{(T,\cdot)}G^{k}(t,x;T,\cdot)\right)(\xi)=\psi(\xi)\hat{G}^{k}(t,x;T,{\xi})-\partial_{T}\hat{G}^{k}(t,x;T,{\xi}),

where $\psi(\xi)$ is the characteristic exponent of the Lévy process with coefficients $\gamma_{0}$ , $a_{0}$ and $\nu_{0}(dz)$ , that is

\psi(\xi)=i\xi(r+\gamma_{0})+a_{0}(-\xi^{2}-i\xi)-\gamma_{0}-\int_{\mathbb{R}}\nu_{0}(dz)(e^{z}-1-z)i\xi+\int_{\mathbb{R}}\nu_{0}(dz)(e^{iz\xi}-1-iz\xi).

(2.21)

Thus the solution (in the Fourier space) to problems (2.14) and (2.16) is given by

\begin{split}\hat{G}^{0}(t,x;T,\xi)&=e^{i\xi x+(T-t)\psi(\xi)},\\ \hat{G}^{k}(t,x;T,\xi)&=-\int_{t}^{T}e^{\psi(\xi)(T-s)}\mathcal{F}\left(\sum_{h=1}^{k}\left(\tilde{L}_{h}^{(s,\cdot)}-\tilde{L}_{h-1}^{(s,\cdot)}\right)G^{k-h}(t,x;s,\cdot)\right)(\xi)ds,\qquad k\geq 1.\end{split}

(2.22)

Now we consider the general framework and in particular we drop the assumption on the existence of the fundamental solution of $L_{0}$ : in this case, we define the $n$ th-order approximation of the characteristic function $\hat{\Gamma}$ as in (2.13), with $\hat{G}^{k}$ given by (2.22). We also notice that

	$\displaystyle\mathcal{F}\left(\left(\tilde{L}_{h}^{(s,\cdot)}-\tilde{L}_{h-1}^{(s,\cdot)}\right)u(s,\cdot)\right)({\xi})=$		(2.23)
	$\displaystyle\qquad\left(a_{h}h(h-1)(-i\partial_{\xi}-\bar{x})^{h-2}+a_{h}(-i\partial_{\xi}-\bar{x})^{h-1}\left(-2hi\xi+(-i\partial_{\xi}-\bar{x})(-\xi^{2}-i\xi)+h\right)\right)\hat{u}(s,\xi)$		(2.24)
	$\displaystyle\qquad-\left(\gamma_{h}h(-i\partial_{\xi}-\bar{x})^{h-1}-\gamma_{h}(-i\partial_{\xi}-\bar{x})^{h}\left(i\xi-1\right)\right)\hat{u}(s,\xi)$		(2.25)
	$\displaystyle\qquad+\int_{\mathbb{R}}\nu_{h}(dz)(e^{z}-1-z)\left(h(-i\partial_{\xi}-\bar{x})^{h-1}-(-i\partial_{\xi}-\bar{x})^{h}i\xi\right)\hat{u}(s,\xi)$
	$\displaystyle\qquad+\int_{\mathbb{R}}\nu_{h}(dz)\left((-i\partial_{y}-z-\bar{x})^{h}e^{i\xi z}-(-i\partial_{y}-\bar{x})^{h}+z\left(h(-i\partial_{\xi}-\bar{x})^{h-1}-(-i\partial_{\xi}-\bar{x})^{h}i\xi\right)\right)\hat{u}(s,\xi).$		(2.26)

Remark 2.2.

In case the coefficients $\gamma$ , $\sigma$ , $\nu$ depend on time, the solutions to the Cauchy problems are similar:

\begin{split}\hat{G}^{0}(t,x;T,\xi)&=e^{i\xi x}e^{\int_{t}^{T}\psi(s,\xi)ds},\\ \hat{G}^{k}(t,x;T,\xi)&=-\int_{t}^{T}e^{\int_{s}^{T}\psi(\tau,\xi)d\tau}\mathcal{F}\left(\sum_{h=1}^{k}\left(\tilde{L}_{h}^{(s,\cdot)}(s)-\tilde{L}_{h-1}^{(s,\cdot)}(s)\right)G^{k-h}(t,x;s,\cdot)\right)(\xi)ds,\end{split}

(2.27)

with

	$\displaystyle\psi(s,\xi)=i\xi(r+\gamma_{0}(s))+a_{0}(s)(-\xi^{2}-i\xi)-\int_{\mathbb{R}}\nu_{0}(s,dz)(e^{z}-1-z)i\xi+\int_{\mathbb{R}}\nu_{0}(s,dz)(e^{iz\xi}-1-iz\xi),$		(2.28)
	$\displaystyle\tilde{L}_{h}^{(s,y)}(s)-\tilde{L}_{h-1}^{(s,y)}(s)=a_{h}(s)h(h-1)(y-\bar{x})^{h-2}+a_{h}(s)(y-\bar{x})^{h-1}\left(2h\partial_{y}+(y-\bar{x})(\partial_{yy}+\partial_{y})+h\right)$		(2.29)
	$\displaystyle\qquad\qquad\qquad\qquad\qquad-\gamma_{h}(s)h(y-\bar{x})^{h-1}-\gamma_{h}(s)(y-\bar{x})^{h}\left(\partial_{y}+1\right)$		(2.30)
	$\displaystyle\qquad\qquad\qquad\qquad\qquad+\int_{\mathbb{R}}\nu_{h}(s,dz)(e^{z}-1-z)\left(h(y-\bar{x})^{h-1}+(y-\bar{x})^{h}\partial_{y}\right)$
	$\displaystyle\qquad\qquad\qquad\qquad\qquad+\int_{\mathbb{R}}\bar{\nu}_{h}(s,dz)\left((y+z-\bar{x})^{h}e^{z\partial_{y}}-(y-\bar{x})^{h}-z\left(h(y-\bar{x})^{h-1}-(y-\bar{x})^{h}\partial_{y}\right)\right).$		(2.31)

From these results one can already see that the dependency on $x$ comes in through $e^{i\xi x}$ and after taking derivatives the dependency on $x$ will take the form $(x-\bar{x})^{m}e^{i\xi x}$ : this fact will be crucial in our analysis.

Example 2.3.

To see the above dependency more explicitly for the second-order approximation of the characteristic function we consider, for ease of notation, a simplified model: a one-dimensional local Lévy model where the log-price solves the SDE

dX_{t}=\mu(X_{t})dt+\sigma(X_{t})dW_{t}+\int_{\mathbb{R}}d\tilde{N}_{t}(dz)z.

(2.32)

This model is a simplification of the original model, since we consider only a local volatility function, and no local default or state-dependent Lévy measure. Thus only a Taylor expansion of the local volatility coefficient is used. However, the dependency that we will see generalizes in the same way to the local default and state-dependent measure. By the martingale condition we have

\mu(x)=r-a(x)-\int_{\mathbb{R}}\nu(dz)(e^{z}-1),

and therefore the Kolmogorov operator of (2.32) reads

	$\displaystyle Lu(t,x)=$	$\displaystyle\ \partial_{t}u(t,x)+r\partial_{x}u(t,x)+a(t,x)(\partial_{xx}-\partial_{x})u(t,x)$		(2.33)
		$\displaystyle-\int_{\mathbb{R}}\nu(dz)(e^{z}-1)+\int_{\mathbb{R}}\nu(dz)\left(u(t,x+z)-u(t,x)\right).$		(2.34)

In this case, we have the following explicit approximation formulas for the characteristic function $\hat{\Gamma}(t,x;T,{\xi})$ :

\hat{\Gamma}(t,x;T,{\xi})\ \approx\ \hat{\Gamma}^{(n)}(t,x;T,{\xi}):=e^{i\xi x+(T-t)\psi(\xi)}\sum_{k=0}^{n}\hat{F}^{k}(t,x;T,\xi),\qquad n\geq 0,

(2.35)

with

\psi(\xi)=ir\xi-a_{0}(\xi^{2}+i\xi)-\int_{\mathbb{R}}\nu(dz)(e^{z}-1)i\xi+\int_{\mathbb{R}}\nu(dz)\left(e^{iz\xi}-1\right),

and

\displaystyle\hat{F}^{k}(t,x;T,\xi)

\displaystyle=\sum_{h=0}^{k}g^{(k)}_{h}(T-t,\xi)(x-\bar{x})^{h};

(2.36)

here, for $k=0,1,2$ , we have

$\displaystyle g^{(0)}_{0}(s,\xi)=$	$\displaystyle\ 1,$	(2.37)
$\displaystyle g^{(1)}_{0}(s,\xi)=$	$\displaystyle\ a_{1}s^{2}(\xi^{2}+i\xi)\frac{i}{2}\psi^{\prime}(\xi),$	(2.38)
$\displaystyle g^{(1)}_{1}(s,\xi)=$	$\displaystyle\,-a_{1}s(\xi^{2}+i\xi),$	(2.39)
$\displaystyle g^{(2)}_{0}(s,\xi)=$	$\displaystyle\ \frac{1}{2}s^{2}a_{2}\xi(i+\xi)\psi^{\prime\prime}(\xi)-\frac{1}{6}s^{3}\xi(i+\xi)(a_{1}^{2}(i+2\xi)\psi^{\prime}(\xi)-2a_{2}\psi^{\prime}(\xi)^{2}+a_{1}^{2}\xi(i+\xi)\psi^{\prime\prime}(\xi))$	(2.40)
	$\displaystyle-\frac{1}{8}s^{4}a_{1}^{2}\xi^{2}(i+\xi)^{2}\psi^{\prime}(\xi)^{2},$	(2.41)
$\displaystyle g^{(2)}_{1}(s,\xi)=$	$\displaystyle\ \frac{1}{2}s^{2}\xi(i+\xi)(a_{1}^{2}(1-2i\xi)+2ia_{2}\psi^{\prime\prime}(\xi))-\frac{1}{2}s^{3}ia_{1}^{2}\xi^{2}(i+\xi)^{2}\psi^{\prime\prime}(\xi),$	(2.42)
$\displaystyle g^{(2)}_{2}(s,\xi)=$	$\displaystyle\,-a_{2}s\xi(i+\xi)+\frac{1}{2}s^{2}a_{1}^{2}\xi^{2}(i+\xi)^{2}.$	(2.43)

Using the notation from above, we can write in the same way the approximation formulas for the general case. Here we present the results for $k=0,1$ , since higher-order formulas are too long to include. For the full formula we refer to Appendix B. We have:

$\displaystyle g_{0}^{(0)}(s,\xi)=$	$\displaystyle 1,$	(2.44)
$\displaystyle g_{0}^{(1)}(s,\xi)=$	$\displaystyle\frac{i}{2}a_{1}s^{2}({\xi}^{2}+i{\xi})\psi^{\prime}({\xi})+\frac{1}{2}\gamma_{1}s^{2}(i+{\xi})\psi^{\prime}({\xi})-\frac{1}{2}\int_{\mathbb{R}}\nu_{1}(dz)(e^{z}-1-z)s^{2}{\xi}\psi^{\prime}({\xi})$	(2.45)
	$\displaystyle-\frac{1}{2}\int_{\mathbb{R}}\nu_{1}(dz)(ie^{i\xi z}-i+{\xi}z)s^{2}\psi^{\prime}({\xi}),$	(2.46)
$\displaystyle g_{1}^{(1)}(s,{\xi})=$	$\displaystyle-a_{1}s(\xi^{2}+i\xi)+\gamma_{1}si(i+{\xi})-\int_{\mathbb{R}}\nu_{1}(dz)(e^{z}-1-z)s{\xi}i$	(2.47)
	$\displaystyle+\int_{\mathbb{R}}\nu_{1}(dz)(e^{i{\xi}z}-1-{\xi}iz)s.$	(2.48)

Remark 2.4.

From (2.35)-(2.36) and (2.49) we clearly see that the approximation of order $n$ is a function of the form

\hat{\Gamma}^{(n)}(t,x;T,{\xi}):=e^{i\xi x}\sum_{k=0}^{n}(x-\bar{x})^{k}g_{n,k}(t,T,{\xi}),

(2.49)

where the coefficients $g_{n,k}$ , with $0\leq k\leq n$ , depend only on $t,T$ and ${\xi}$ , but not on $x$ . The approximation formula can thus always be split into a sum of products of functions depending only on $\xi$ and functions that are linear combinations of $(x-\bar{x})^{m}e^{i\xi x}$ , $m\in\mathbb{N}_{0}$ .

3 Bermudan option valuation

A Bermudan option is a financial contract in which the holder can exercise at a predetermined finite set of exercise moments prior to maturity, and the holder of the option receives a payoff when exercising. Consider a Bermudan option with a set of $M$ exercise moments $\{t_{1},...,t_{M}\}$ , with $0\leq t_{1}<t_{2}<\cdots<t_{M}=T$ . When the option is exercised at time $t_{m}$ the holder receives the payoff $\Phi\left(t_{m},S_{t_{m}}\right)$ . Recalling (2.3), the no-arbitrage value of the Bermudan option at time $t$ is

\displaystyle v\left(t,X_{t}\right)={\mathds{1}}_{\{\zeta>t\}}\sup_{{\tau}\in\mathcal{T}_{t}}E\left[e^{-\int_{t}^{{\tau}}\left(r+\gamma(s,X_{s})\right)ds}{\varphi}({\tau},X_{{\tau}})|X_{t}\right],

(3.50)

where ${\varphi}(t,x)=\Phi(t,e^{x})$ and $\mathcal{T}_{t}$ is the set of all $\mathcal{G}$ -stopping times taking values in $\{t_{1},...,t_{M}\}\cap[t,T]$ . For a Bermudan Put option with strike price $K$ , we simply have ${\varphi}(t,x)=\left(K-e^{x}\right)^{+}$ . By the dynamic programming approach, the option value can be expressed by a backward recursion as

v(t_{M},x)={\mathds{1}}_{\{\zeta>t_{M}\}}{\varphi}(t_{M},x)

and

\begin{cases}c(t,x)=E\left[e^{\int_{t}^{t_{m}}\left(r+\gamma(s,X_{s})\right)ds}v(t_{m},X_{t_{m}})|X_{t}=x\right],\qquad&t\in[t_{m-1},t_{m}[\\ v(t_{m-1},x)={\mathds{1}}_{\{\zeta>t_{m-1}\}}\max\{{\varphi}(t_{m-1},x),c(t_{m-1},x)\},\qquad&m\in\{2,\dots,M\}.\end{cases}

(3.51)

In the above notation $v(t,x)$ is the option value and $c(t,x)$ is the so-called continuation value. The option value is set to be $v(t,x)=c(t,x)$ for $t\in\,]t_{m-1},t_{m}[$ , and, if $t_{1}>0$ , also for $t\in[0,t_{1}[$ .

Remark 3.5.

Since the payoff of a Call option grows exponentially with the log-stock price, this may introduce significant cancellation errors for large domain sizes. For this reason we price Put options only using our approach and we employ the well-known Put-Call parity to price Calls via Puts. This is a rather standard argument (see, for instance, [17]).

3.1 An algorithm for pricing Bermudan Put options

The COS method proposed by [5] is based on the insight that the Fourier-cosine series coefficients of $\Gamma(t,x;T,dy)$ (and therefore also of option prices) are closely related to the characteristic function of the underlying process, namely the following relationship holds:

\int_{a}^{b}e^{i\frac{k\pi}{b-a}}\Gamma(t,x;T,dy)\approx\hat{\Gamma}\left(t,x;T,\frac{k\pi}{b-a}\right).

The COS method provides a way to calculating expected values (integrals) of the form

v(t,x)=\int_{\mathbb{R}}{\varphi}(T,y)\Gamma(t,x;T,dy),

and it consists of three approximation steps:

1.

In the first step we truncate the infinite integration range to $[a,b]$ to obtain approximation $v_{1}$ :

$v_{1}(t,x):=\int_{a}^{b}{\varphi}(T,y)\Gamma(t,x;T,dy).$

We assume this can be done due to the rapid decay of the distribution at infinity.

In the second step we replace the distribution with its cosine expansion and we get

v_{1}(t,x):=\frac{b-a}{2}\sideset{}{{}^{\prime}}{\sum}_{k=0}^{\infty}A_{k}(t,x;T)V_{k}(T),

where $\sideset{}{{}^{\prime}}{\sum}$ indicates that the first term in the summation is weighted by one-half and

	$\displaystyle A_{k}(t,x;T)$	$\displaystyle=\frac{2}{b-a}\int_{a}^{b}\cos\left(k\pi\frac{y-a}{b-a}\right)\Gamma(t,x;T,dy),$		(3.52)
	$\displaystyle V_{k}(T)$	$\displaystyle=\frac{2}{b-a}\int_{a}^{b}\cos\left(k\pi\frac{y-a}{b-a}\right){\varphi}(T,y)dy,$		(3.53)

are the Fourier-cosine series coefficients of the distribution and of the payoff function at time $T$ respectively. Due to the rapid decay of the Fourier-cosine series coefficients, we truncate the series summation and obtain approximation $v_{2}$ :

v_{2}(t,x):=\frac{b-a}{2}\sideset{}{{}^{\prime}}{\sum}_{k=0}^{N-1}A_{k}(t,x;T)V_{k}(T).

In the third step we use the fact that the coefficients $A_{k}$ can be rewritten using the truncated characteristic function:

A_{k}(t,x;T)=\frac{2}{b-a}\textnormal{Re}\left(e^{-ik\pi\frac{a}{b-a}}\int_{a}^{b}e^{i\frac{k\pi}{b-a}y}\Gamma(t,x;T,dy)\right).

The finite integration range can be approximated as

\int_{a}^{b}e^{i\frac{k\pi}{b-a}y}\Gamma(t,x;T,dy)\approx\int_{\mathbb{R}}e^{i\frac{k\pi}{b-a}y}\Gamma(t,x;T,dy)=\hat{\Gamma}\left(t,x;T,\frac{k\pi}{b-a}\right).

Thus in the last step we replace $A_{k}$ by its approximation:

\frac{2}{b-a}\textnormal{Re}\left(e^{-ik\pi\frac{a}{b-a}}\hat{\Gamma}\left(t,x;T,\frac{k\pi}{b-a}\right)\right),

(3.54)

and obtain approximation $v_{3}$ :

v_{3}(t,x):=\sideset{}{{}^{\prime}}{\sum}_{k=0}^{N-1}\textnormal{Re}\left(e^{-ik\pi\frac{a}{b-a}}\hat{\Gamma}\left(t,x;T,\frac{k\pi}{b-a}\right)\right)V_{k}(T).

(3.55)

Next we go back to the Bermudan Put pricing problem. Remembering that the expected value $c(t,x)$ in (3.51) can be rewritten in integral form as in (2.7), we have

\displaystyle c(t,x)=e^{-r(t_{m}-t)}\int_{\mathbb{R}}v(t_{m},y)\Gamma(t,x;t_{m},dy),\qquad t\in[t_{m-1},t_{m}[.

(3.56)

Then we use the Fourier-cosine expansion (3.55), so that we get the approximation:

		$\displaystyle\hat{c}(t,x)=e^{-r(t_{m}-t)}\sideset{}{{}^{\prime}}{\sum}_{k=0}^{N-1}\textnormal{Re}\left(e^{-ik\pi\frac{a}{b-a}}\hat{\Gamma}\left(t,x;t_{m},\frac{k\pi}{b-a}\right)\right)V_{k}(t_{m}),\qquad t\in[t_{m-1},t_{m}[$		(3.57)
		$\displaystyle V_{k}(t_{m})=\frac{2}{b-a}\int_{a}^{b}\cos\left(k\pi\frac{y-a}{b-a}\right)\max\{{\varphi}(t_{m},y),c(t_{m},y)\}dy,$		(3.58)

with ${\varphi}(t,x)=\left(K-e^{x}\right)^{+}$ .

Next we recover the coefficients $\left(V_{k}(t_{m})\right)_{k=0,1,...,N-1}$ from $\left(V_{k}(t_{m+1})\right)_{k=0,1,...,N-1}$ . To this end, we split the integral in the definition of $V_{k}(t_{m})$ into two parts using the early-exercise point $x_{m}^{*}$ , which is the point where the continuation value is equal to the payoff, i.e. $c(t_{m},x_{m}^{*})={\varphi}(t_{m},x_{m}^{*})$ ; thus we have

V_{k}(t_{m})=F_{k}(t_{m},x_{m}^{*})+C_{k}(t_{m},x_{m}^{*}),\qquad m=M-1,M-2,...,1,

where

\begin{split}F_{k}(t_{m},x_{m}^{*})&:=\frac{2}{b-a}\int_{a}^{x_{m}^{*}}{\varphi}(t_{m},y)\cos\left(k\pi\frac{y-a}{b-a}\right)dy,\\ C_{k}(t_{m},x_{m}^{*})&:=\frac{2}{b-a}\int_{x_{m}^{*}}^{b}c(t_{m},y)\cos\left(k\pi\frac{y-a}{b-a}\right)dy,\end{split}

(3.59)

and $V_{k}(t_{M})=F_{k}(t_{M},\log K).$

Remark 3.6.

Since we have a semi-analytic formula for $\hat{c}(t_{m},x)$ , we can easily find the derivatives with respect to $x$ and use Newton’s method to find the point $x_{m}^{*}$ such that $c(t_{m},x_{m}^{*})={\varphi}(t_{m},x_{m}^{*})$ . A good starting point for the Newton method is $\log K$ , since $x_{m}^{*}\leq\log K$ .

The coefficients $F_{k}(t_{m},x_{m}^{*})$ can be computed analytically using $x_{m}^{*}\leq\log K$ , so that we have

	$\displaystyle F_{k}(t_{m},x_{m}^{*})$	$\displaystyle=\frac{2}{b-a}\int_{a}^{x_{m}^{*}}(K-e^{y})\cos\left(k\pi\frac{y-a}{b-a}\right)dy$		(3.60)
		$\displaystyle=\frac{2}{b-a}K\Psi_{k}(a,x_{m}^{})-\frac{2}{b-a}\chi_{k}(a,x_{m}^{}),$		(3.61)

where

$\displaystyle\chi_{k}(a,x_{m}^{*})$	$\displaystyle=\int_{a}^{x_{m}^{*}}e^{y}\cos\left(k\pi\frac{y-a}{b-a}\right)dy$	(3.62)
	$\displaystyle=\frac{1}{1+\left(\frac{k\pi}{b-a}\right)^{2}}\left(e^{x_{m}^{}}\cos\left(k\pi\frac{x_{m}^{}-a}{b-a}\right)-e^{a}+\frac{k\pi e^{x_{m}^{}}}{b-a}\sin\left(k\pi\frac{x_{m}^{}-a}{b-a}\right)\right),$	(3.63)
$\displaystyle\Psi_{k}(a,x_{m}^{*})$	$\displaystyle=\int_{a}^{x_{m}^{}}\cos\left(k\pi\frac{y-a}{b-a}\right)dy=\begin{cases}\frac{b-a}{k\pi}\sin\left(k\pi\frac{x_{m}^{}-a}{b-a}\right),\qquad&k\neq 0,\\ x_{m}^{*}-a,&k=0.\end{cases}$	(3.64)

On the other hand, by inserting the approximation (3.57) for the continuation value into the formula for $C_{k}(t_{m},x_{m}^{*})$ have the following coefficients $\hat{C}_{k}$ for $m=M-1,M-2,...,1$ :

\hat{C}_{k}(t_{m},x_{m}^{*})=\frac{2e^{-r(t_{m+1}-t_{m})}}{b-a}\sideset{}{{}^{\prime}}{\sum}_{j=0}^{N-1}V_{j}(t_{m+1})\int_{x_{m}^{*}}^{b}\mathrm{Re}\left(e^{-ij\pi\frac{a}{b-a}}\hat{\Gamma}\left(t_{m},x;t_{m+1},\frac{j\pi}{b-a}\right)\right)\cos\left(k\pi\frac{x-a}{b-a}\right)dx.

(3.65)

Thus the algorithm for pricing Bermudan options can then be summarized as follows:

Figure 1: Algorithm 3.1: Bermudan option valuation

1. For

k=0,1,...,N-1

: • At time

t_{M}

, the coefficients are exact:

V_{k}(t_{M})=F_{k}(t_{M},\log K)

, as in (3.59). 2. For

m=M-1

to 1: • Determine the early-exercise point

x_{m}^{*}

using Newton’s method; • Compute

\hat{V}_{k}(t_{m})

using formula

\hat{V}_{k}(t_{m}):=F_{k}(t_{m},x_{m}^{*})+\hat{C}_{k}(t_{m},x_{m}^{*})

, (3.59) and (3.65). Use an FFT for the continuation value (see Section 3.2). 3. Final step: using

\hat{V}_{k}(t_{1})

determine the option price

\hat{v}(0,x)=\hat{c}(0,x)

using (3.57).

3.2 An efficient algorithm for the continuation value

In this section we derive an efficient algorithm for calculating $\hat{C}_{k}(t_{m},x_{m}^{*})$ in (3.65). When considering an exponential Lévy process with constant coefficients as done in [5], the continuation value can be calculated using a Fast Fourier Transform (FFT). This can be done due to the fact that the characteristic function $\hat{\Gamma}(t,x;T,{\xi})$ can be split into a product of a function depending only on $\xi$ and a function of the form $e^{i\xi x}$ . Note that we typically have $\xi=\frac{j\pi}{b-a}$ . The integration over $x$ results in a sum of a Hankel and Toeplitz matrix (with indices $(j+k)$ and $(j-k)$ respectively). The matrix-vector product, with these special matrices, can be transformed into a circular convolution which can be computed using FFTs.

From (2.49) we know that the $n$ th-order approximation of the characteristic function is of the form:

\hat{\Gamma}^{(n)}(t_{m},x;t_{m+1},{\xi})=e^{i\xi x}\sum_{k=0}^{n}(x-\bar{x})^{k}g_{n,k}(t_{m},t_{m+1},{\xi}),

where the coefficients $g_{n,k}(t,T,{\xi})$ , with $0\leq k\leq n$ , depend only on $t,T$ and ${\xi}$ , but not on $x$ . Using (2.49) we write the continuation value as:

\displaystyle\hat{C}_{k}(t_{m},x_{m}^{*})

\displaystyle=\sum_{h=0}^{n}e^{-r(t_{m+1}-t_{m})}\sideset{}{{}^{\prime}}{\sum}_{j=0}^{N-1}\mathrm{Re}\left(V_{j}(t_{m})g_{n,h}\left(t_{m},t_{m+1},\frac{j\pi}{b-a}\right)M^{h}_{k,j}(x_{m}^{*},b)\right),

(3.66)

where we have interchanged the sums and integral and defined:

\displaystyle M_{k,j}^{h}(x_{m}^{*},b)=\frac{2}{b-a}\int_{x_{m}^{*}}^{b}e^{ij\pi\frac{x-a}{b-a}}(x-\bar{x})^{h}\cos\left(k\pi\frac{x-a}{b-a}\right)dx

(3.67)

This can be written in vectorized form as:

\displaystyle\mathbb{\hat{C}}_{k}(t_{m},x_{m}^{*})=\sum_{h=1}^{n}e^{-r(t_{m+1}-t_{m})}\mathrm{Re}\left(\mathbb{V}(t_{m+1})\mathcal{M}^{h}(x_{m}^{*},b)\Lambda^{h}\right),

(3.68)

where $\mathbb{V}(t_{m+1})$ is the vector $[V_{0}(t_{m+1}),...,V_{N-1}(t_{m+1})]^{T}$ and $\mathcal{M}^{h}(x_{m}^{*},b)\Lambda^{h}$ is a matrix-matrix product with $\mathcal{M}^{h}$ being a matrix with elements $\{M_{k,j}^{h}\}_{k,j=0}^{N-1}$ and $\Lambda^{h}$ is a diagonal matrix with elements

g_{n,h}\Big{(}t_{m},t_{m+1},\frac{j\pi}{b-a}\Big{)},\qquad j=0,\dots,N-1.

We have the following theorem for calculating a generalized form of the integral in (3.67) which is used in the calculation of the continuation value.

Theorem 3.7.

The matrix $\mathcal{M}$ with elements $\{M_{k,j}\}_{k,j=0}^{N-1}$ such that:

\displaystyle M_{k,j}=\int e^{jx}\cos(kx)x^{m}dx,

(3.69)

consists of sums of Hankel and Toeplitz matrices.

Proof.

Using standard trigonometric identities we can rewrite the integral as:

	$\displaystyle M_{k,j}$	$\displaystyle=\int\cos(jx)\cos(kx)x^{m}dx+i\int\sin(jx)\cos(kx)x^{m}dx$		(3.70)
		$\displaystyle=M_{k,j}^{H}+iM_{k,j}^{T},$		(3.71)

where we have defined:

	$\displaystyle M_{k,j}^{H}$	$\displaystyle=\frac{1}{2}\int\cos((j+k)x)x^{m}dx+\frac{1}{2}\int\sin((j+k)x)x^{m}dx,$		(3.72)
	$\displaystyle M_{k,j}^{T}$	$\displaystyle=\frac{1}{2}\int\cos((j-k)x)x^{m}dx+\frac{1}{2}\int\sin((j-k)x)x^{m}dx.$		(3.73)

The following holds:

$\displaystyle\int\cos(nx)x^{m}dx=$	$\displaystyle\,\frac{1}{n}x^{m}\sin(nx)+\sum_{i=1}^{\lceil m/2\rceil}(-1)^{i+1}\frac{\prod_{j=0}^{2i-2}(m-j)}{n^{2i}}\cos(nx)x^{m-(2i-1)}$	(3.74)
	$\displaystyle-\sum_{i=1}^{\lfloor m/2\rfloor}(-1)^{i+1}\frac{\prod_{j=0}^{2i-1}(m-j)}{n^{2i+1}}\sin(nx)x^{m-2i},$	(3.75)
$\displaystyle\int\sin(nx)x^{m}dx=$	$\displaystyle\,-\frac{1}{n}x^{m}\cos(nx)+\sum_{i=1}^{\lceil m/2\rceil}(-1)^{i+1}\frac{\prod_{j=0}^{2i-2}(m-j)}{n^{2i}}\sin(nx)x^{m-(2i-1)}$	(3.76)
	$\displaystyle-\sum_{i=1}^{\lfloor m/2\rfloor}(-1)^{i+1}\frac{\prod_{j=0}^{2i-1}(m-j)}{n^{2i+1}}\cos(nx)x^{m-2i}.$	(3.77)

It follows that $\displaystyle\{M_{k,j}^{H}\}_{k,j=0}^{N-1}$ is a Hankel matrix with coefficient $(j+k)$ and $\displaystyle\{M_{k,j}^{T}\}_{k,j=0}^{N-1}$ is a Toeplitz matrix with coefficient $(j-k)$ :

\mathcal{M}_{H}=\begin{pmatrix}M_{0}&M_{1}&M_{2}&\dots&M_{N-1}\\ M_{1}&M_{2}&\dots&&M_{N}\\ \vdots&&&&\vdots\\ M_{N-2}&M_{N-1}&\dots&&M_{2N-3}\\ M_{N-1}&\dots&&M_{2N-3}&M_{2N-2}\end{pmatrix},

\mathcal{M}_{T}=\begin{pmatrix}M_{0}&M_{1}&\dots&M_{N-2}&M_{N-1}\\ M_{-1}&M_{0}&M_{1}&\dots&M_{N-2}\\ \vdots&&\ddots&&\vdots\\ M_{2-N}&\dots&M_{-1}&M_{0}&M_{1}\\ M_{1-N}&M_{2-N}&&M_{-1}&M_{0}\end{pmatrix},

where we have defined

\displaystyle M_{j}=\frac{1}{2}\int\cos(jx)x^{m}dx+\frac{1}{2}\int\sin(jx)x^{m}dx.

(3.78)

∎

From Theorem 3.7 we see that $\mathcal{M}^{h}(x_{m}^{*},b)$ with elements $M_{k,j}^{h}$ consists of a sum of a Hankel and Toeplitz matrix.

Example 3.8.

We derive explicitly the Hankel and Toeplitz matrices for $m=0$ and $m=1$ . We calculate the indefinite integral

\displaystyle M_{k,j}=\frac{2}{b-a}\int e^{ij\pi\frac{x-a}{b-a}}\cos\left(k\pi\frac{x-a}{b-a}\right)(x-\bar{x})^{m}dx.

(3.79)

Suppose $m=0$ , in this case we have $M_{k,j}=M_{k,j}^{H}+M_{k,j}^{T}$ , with:

	$\displaystyle M_{k,j}^{H}$	$\displaystyle=-\frac{i\exp\left(i\frac{(j+k)\pi(x-a)}{b-a}\right)}{\pi(j+k)},$		(3.80)
	$\displaystyle M_{k,j}^{T}$	$\displaystyle=-\frac{i\exp\left(i\frac{(j-k)\pi(x-a)}{b-a}\right)}{\pi(j-k)},$		(3.81)

where $\displaystyle\{M_{k,j}^{H}\}_{k,j=0}^{N-1}$ is a Hankel matrix and $\displaystyle\{M_{k,j}^{T}\}_{k,j=0}^{N-1}$ is a Toeplitz matrix with

\displaystyle M_{j}=\begin{cases}\frac{x}{b-a},\qquad&j=0,\\ \frac{i\exp\left(i\frac{j\pi(x-a)}{b-a}\right)}{\pi j},\quad&j\neq 0.\end{cases}

(3.82)

Suppose $m=1$ , in this case we have:

	$\displaystyle M_{k,j}^{H}$	$\displaystyle=-\frac{a-b}{(j-k)^{2}\pi^{2}}\exp\left(i(j-k)\pi\frac{(x-a)}{b-a}\right)-\frac{x-\bar{x}}{(j-k)\pi}i\exp\left(i(j-k)\pi\frac{(x-a)}{b-a}\right),$		(3.83)
	$\displaystyle M_{k,j}^{T}$	$\displaystyle=-\frac{a-b}{(j+k)^{2}\pi^{2}}\exp\left(i(j+k)\pi\frac{(x-a)}{b-a}\right)-\frac{x-\bar{x}}{(j+k)\pi}i\exp\left(i(j+k)\pi\frac{(x-a)}{b-a}\right),$		(3.84)

where $\displaystyle\{M_{k,j}^{H}\}_{k,j=0}^{N-1}$ is a Hankel matrix and $\displaystyle\{M_{k,j}^{T}\}_{k,j=0}^{N-1}$ is a Toeplitz matrix, with

\displaystyle M_{j}=\begin{cases}\frac{x(x-\bar{x})}{b-a},\qquad&j=0,\\ -\frac{a-b}{j^{2}\pi^{2}}\exp\left(ij\pi\frac{(x-a)}{b-a}\right)-\frac{x-\bar{x}}{j\pi}i\exp\left(ij\pi\frac{(x-a)}{b-a}\right),\quad&j\neq 0.\end{cases}

(3.85)

Remark 3.9.

If we take $\bar{x}=x$ , which is most common in practice, the formulas are simplified significantly and only the case of $m=0$ is relevant. In this case the characteristic function is simply $e^{i\xi x}$ times a sum of terms depending only on $t_{m}$ , $t_{m+1}$ and $\xi=\frac{j\pi}{b-a}$ :

\hat{\Gamma}^{(n)}(t_{m},x;t_{m+1},{\xi})=e^{i\xi x}g_{n,0}(t_{m},t_{m+1},{\xi}).

Using the split into sums of Hankel and Toeplitz matrices we can write the continuation value in matrix form as:

\displaystyle\boldsymbol{\hat{C}}(t_{m},x_{m}^{*})=\sum_{h=0}^{n}e^{-r(t_{m+1}-t_{m})}\mathrm{Re}\left((\mathcal{M}^{h}_{H}+\mathcal{M}^{h}_{T})\boldsymbol{u}^{l}\right),

(3.86)

where $\mathcal{M}^{h}_{H}=\{M^{H,h}_{k,j}(x_{m}^{*},b)\}_{k,j=0}^{N-1}$ is a Hankel matrix and $\mathcal{M}^{l}_{T}=\{M^{T,h}_{k,j}(x_{m}^{*},b)\}_{k,j=0}^{N-1}$ is a Toeplitz matrix and $\boldsymbol{u}^{h}=\{u_{j}^{h}\}_{j=0}^{N-1}$ , with $u_{j}^{h}=g_{n,h}\left(t_{m},t_{m+1},\frac{j\pi}{b-a}\right)V_{j}(t_{m+1})$ and $u_{0}^{h}=\frac{1}{2}g_{n,h}\left(t_{m},t_{m+1},0\right)V_{0}(t_{m+1})$ .

We recall that the circular convolution, denoted by $\circledast$ , of two vectors is equal to the inverse discrete Fourier transform $(\mathcal{D}^{-1})$ of the products of the forward DFTs, $\mathcal{D}$ , i.e.:

\mathbb{x}\circledast\mathbb{y}=\mathcal{D}^{-1}\{\mathcal{D}(\mathbb{x})\cdot\mathcal{D}(\mathbb{y})\}.

For Hankel and Toeplitz matrices we have the following result:

Theorem 3.10.

For a Toeplitz matrix $\mathcal{M}_{T}$ , the product $\mathcal{M}_{T}\mathbb{u}$ is equal to the first $N$ elements of $\mathbb{m}_{T}\circledast\mathbb{u}_{T}$ , where $\mathbb{m}_{T}$ and $\mathbb{u}_{T}$ are $2N$ vectors defined by

	$\displaystyle\mathbb{m}_{T}=[M_{0},M_{-1},M_{-2},...,M_{1-N},0,M_{N-1},M_{N-2},...,M_{1}]^{T},$		(3.87)
	$\displaystyle\mathbb{u}_{T}=[u_{0},u_{1},...,u_{N-1},0,...,0]^{T}.$		(3.88)

For a Hankel matrix $\mathcal{M}_{H}$ , the product $\mathcal{M}_{H}\mathbb{u}$ is equal to the first $N$ elements of $\mathbb{m_{H}}\circledast\mathbb{u_{H}}$ in reversed order, where $\mathbb{m}_{H}$ and $\mathbb{u}_{H}$ are $2N$ vectors defined by

	$\displaystyle\mathbb{m}_{H}=[M_{2N-1},M_{2N-2},...,M_{1},M_{0}]^{T}$		(3.89)
	$\displaystyle\mathbb{u}_{H}=[0,...,0,u_{0},u_{1},...,u_{N-1}]^{T}.$		(3.90)

Summarizing, we can calculate the continuation value $\mathbb{\hat{C}}(t_{m},x_{m}^{*})$ using the algorithm in Figure 2.

Figure 2: Algorithm 3.2: Computation of

\mathbb{\hat{C}}(t_{m},x_{m}^{*})

1. For

h=0,...,n

: • Compute

M^{h}_{j}(x_{1},x_{2})

• Construct

\mathbb{m}^{h}_{H}

and

\mathbb{m}^{h}_{T}

• Compute

\boldsymbol{u}^{h}(t_{m})=\{u_{j}^{h}\}_{j=0}^{N-1}

• Construct

\mathbb{u}^{h}_{T}

by padding

N

zeros to

\boldsymbol{u}^{h}(t_{m})

•

\mathbb{MTu}^{h}=

the first

N

elements of

\mathcal{D}^{-1}\{\mathcal{D}(\mathbb{m}_{T}^{h})\cdot\mathcal{D}(\mathbb{u}_{T}^{h})\}

•

\mathbb{MHu}^{h}=

reverse

\{

the first

N

elements of

\mathcal{D}^{-1}\{\mathcal{D}(\mathbb{m}_{H}^{h})\cdot\mathbb{sgn}\cdot\mathcal{D}(\mathbb{u}_{T}^{h})\}\}

2. Compute the continuation value using

\mathbb{\hat{C}}(t_{m},x_{m}^{*})=\sum\limits_{h=0}^{n}e^{-r(t_{m+1}-t_{m})}\mathrm{Re}(\mathbb{MTu}^{h}+\mathbb{MHu}^{h})

The continuation value requires five DFTs for each $h=0,...,n$ , and a DFT is calculated using the FFT. In practice it is most common to have $\bar{x}=x$ and in this case we only need five FFTs. The computation of $F_{k}(t_{m},x_{m}^{*})$ is linear in $N$ . The overall complexity of the method is dominated by the computation of $\hat{C}(t_{m},x_{m}^{*})$ , whose complexity is $O(N\log_{2}N)$ with the FFT. The complexity of the calculation for option value at time 0 is $O(N)$ . If we have a Bermudan option with $M$ exercise dates, the overall complexity will be $O((M-1)N\log_{2}N)$ .

Remark 3.11 (American options).

The prices of American options can be obtained by applying a Richardson extrapolation (see, for instance, [9]) on the prices of a few Bermudan options with a small number of exercise dates. Let $v_{M}$ denote the value of a Bermudan option with maturity $T$ and a number $M$ of early exercise dates that are $\tfrac{T}{M}$ years apart. Then, for any $d\in\mathbb{N}$ , the following 4-point Richardson extrapolation scheme

\frac{1}{21}\left(64v_{2^{d+3}}-56v_{2^{d+2}}+14v_{2^{d+1}}-v_{2^{d}}\right)

gives an approximation of the corresponding American option price.

Remark 3.12 (The Greeks).

The approximation method can also be used to calculate the Greeks at almost no additional cost. In the case of $\bar{x}=x$ , we have the following approximation formulas for Delta and Gamma:

$\displaystyle\hat{\Delta}=$	$\displaystyle\ e^{-r(t_{1}-t_{0})}\sideset{}{{}^{\prime}}{\sum}_{k=0}^{N-1}\textnormal{Re}\left(e^{ik\pi\frac{x-a}{b-a}}\left(\frac{ik\pi}{b-a}g_{n,0}\left(t_{0},t_{1},\frac{k\pi}{b-a}\right)+g_{n,1}\left(t_{0},t_{1},\frac{k\pi}{b-a}\right)\right)\right)\hat{V}_{k}(t_{1}),$	(3.91)
$\displaystyle\hat{\Gamma}=$	$\displaystyle\ e^{-r(t_{1}-t_{0})}\sideset{}{{}^{\prime}}{\sum}_{k=0}^{N-1}\textnormal{Re}\bigg{(}e^{ik\pi\frac{x-a}{b-a}}\bigg{(}-\frac{ik\pi}{b-a}g_{n,0}\left(t_{0},t_{1},\frac{k\pi}{b-a}\right)-g_{n,1}\left(t_{0},t_{1},\frac{k\pi}{b-a}\right)$	(3.92)
	$\displaystyle+2\frac{ik\pi}{b-a}g_{n,1}\left(t_{0},t_{1},\frac{k\pi}{b-a}\right)+\left(\frac{ik\pi}{b-a}\right)^{2}g_{n,0}\left(t_{0},t_{1},\frac{k\pi}{b-a}\right)+2g_{n,2}\left(t_{0},t_{1},\frac{k\pi}{b-a}\right)\bigg{)}\bigg{)}\hat{V}_{k}(t_{1}).$	(3.93)

4 Error estimates

The error in our approximation consists of the error of the COS method and the error in the adjoint expansion of the characteristic function. The error of the COS method depends on the truncation of the integration range $[a,b]$ and the truncation of the infinite summation of the Fourier-cosine expansion by $N$ . The density rapidly decays to zero as $y\rightarrow\pm\infty$ . Then the overall error can be bounded as follows:

\epsilon_{1}(x;N,[a,b])\leq Q\left|\int_{\mathbb{R}\backslash[a,b]}\Gamma(t,x;T,dy)\right|+\left|\frac{P}{(N-1)^{\beta-1}}\right|,

where $P$ and $Q$ are constants not depending on $N$ or $[a,b]$ and $\beta\geq n\geq 1$ , with $n$ being the algebraic index of convergence of the cosine series coefficients. For a sufficiently large integration interval $[a,b]$ , the overall error is dominated by the series truncation error, which converges exponentially. The error in the backward propagation of the coefficients $V_{k}(t_{m})$ is defined as $\epsilon_{2}(k,t_{m}):=V_{k}(t_{m})-\hat{V}_{k}(t_{m})$ . With $[a,b]$ sufficiently large and a probability density function in $C^{\infty}([a,b])$ , the error $\epsilon_{1}(k,t_{m})$ converges exponentially in $N$ . For a detailed derivation on the error of the COS method see [4] and [5].

We now present the error estimates for the adjoint expansion of the characteristic function at orders zero and one. We consider for simplicity a model with time-independent coefficients

X_{t}=x+\int_{0}^{t}\mu(X_{s})ds+\int_{0}^{t}\sigma(X_{s})dW_{s}+\int_{0}^{t}\int_{\mathbb{R}}\eta(X_{s-})zd\tilde{N}(s,dz),

(4.94)

where we have defined as usual $d\tilde{N}(t,dz)=dN(t,dz)-\nu(dz)dt$ . This model is similar to the model we considered initially in (2.1); only now we deal with slightly simplified version and assume that the dependency on $X_{t}$ in the measure can be factored out, which is often enough the case.

Let $\tilde{X}_{t}$ be the 0th-order approximation of the model in (4.94) with $\bar{x}=x$ , that is

\displaystyle\tilde{X}_{t}=x+\int_{0}^{t}\mu(x)ds+\int_{0}^{t}\sigma(x)dW_{s}+\int_{0}^{t}\int_{\mathbb{R}}\eta(x)zd\tilde{N}(s,dz).

(4.95)

The characteristic exponent of $\tilde{X}_{t}-x$ is

\psi(\xi)=i\xi\mu(x)-\frac{\sigma(x)^{2}}{2}{\xi}^{2}-\eta(x)\int_{\mathbb{R}}\nu(dz)(e^{z}-1-z)i\xi+\eta(x)\int_{\mathbb{R}}\nu(dz)(e^{iz\xi}-1-iz\xi).

(4.96)

Theorem 4.13.

Let $n=0,1$ and assume that the coefficients $\mu,{\sigma},{\eta}$ are continuously differentiable with bounded derivatives up to order $n$ . Let $\hat{\Gamma}^{(n)}(0,x;t,{\xi})$ in (2.13) be the $n$ th-order approximation of the characteristic function. Then, for any $T>0$ there exists a positive constant $C$ that depends only on $T$ , on the norms of the coefficients and on the Lévy measure $\nu$ , such that

\displaystyle\left|\hat{\Gamma}(0,x;t,{\xi})-\hat{\Gamma}^{(n)}(0,x;t,{\xi})\right|

\displaystyle\leq C\left(1+|{\xi}|^{1+3n}\right)t^{n+1},\qquad t\in[0,T],\ {\xi}\in{\mathbb{R}}.

(4.97)

Proof.

For the proof we refer to Appendix A. ∎

Remark 4.14.

The proof of Theorem 4.13 can be generalized to obtain error bounds for any $n\in\mathbb{N}$ : however, one can see that, for $n\geq 2$ , the order of convergence improves only in the diffusive part, according to the results proved in [10].

5 Numerical tests

For the numerical examples we use the second-order approximation of the characteristic function. We have found this to be sufficiently accurate by numerical experiments and theoretical error estimates. The formulas for the second-order approximation are simple, making the method easy to implement. For the COS method, unless otherwise mentioned, we use $N=200$ and $L=10$ , where $L$ is the parameter used to define the truncation range $[a,b]$ as follows:

[a,b]:=\left[c_{1}-L\sqrt{c_{2}+\sqrt{c_{4}}},c_{1}+L\sqrt{c_{2}+\sqrt{c_{4}}}\right],

where $c_{n}$ is the $n$ th cumulant of log-price process $\log S$ , as proposed in [4]. The cumulants are calculated using the 0th-order approximation of the characteristic function. A larger $N$ and $L$ has little effect on the price, since a fast convergence is achieved already for small $N$ and $L$ . We compare the approximated values to a 95% confidence interval computed with a Longstaff-Schwartz method with $10^{5}$ simulations and $250$ time steps per year. Furthermore, in the expansion we always use $\bar{x}=x$ .

5.1 Tests under CEV-Merton dynamics

Consider a process under the CEV-Merton dynamics:

dX_{t}=\left(r-a(x)-\lambda\left(e^{m+\delta^{2}/2}-1\right)\right)dt+\sqrt{2a(x)}dW_{t}+\int_{\mathbb{R}}d\tilde{N}_{t}(t,dz)z,

with

	$\displaystyle a(x)=\frac{\sigma_{0}^{2}e^{2(\beta-1)x}}{2},$		(5.98)
	$\displaystyle\nu(dz)=\lambda\frac{1}{\sqrt{2\pi\delta^{2}}}\exp\left(\frac{-(z-m)^{2}}{2\delta^{2}}\right)dz,$		(5.99)
	$\displaystyle\psi(\xi)=-a_{0}(\xi^{2}+i\xi)+ir\xi-i\lambda\left(e^{m+\delta^{2}/2}-1\right)\xi+\lambda\left(e^{mi\xi-\delta^{2}\xi^{2}/2}-1\right).$		(5.100)

We use the following parameters $S_{0}=1$ , $r=5\%$ , $\sigma_{0}=20\%$ , $\beta=0.5$ , $\lambda=30\%$ , $m=-10\%$ , $\delta=40\%$ and compute the European and Bermudan option values.

Table 1: Prices for a European and a Bermudan Put option (expiry

T=0.25

with 3 exercise dates, expiry

T=1

with 10 exercise dates and expiry

T=2

with 20 exercise dates) in the CEV-Merton model for the 2nd-order approximation of the characteristic function, and a Monte Carlo method.

		European		Bermudan
T	K	MC 95% c.i.	Value	MC 95% c.i.	Value
0.25	0.6	0.001240-0.001433	0.001326	0.001243-0.001431	0.001307
	0.8	0.005218-0.005679	0.005493	0.005314-0.005774	0.005421
	1	0.04222-0.04321	0.04275	0.04274-0.04371	0.04304
	1.2	0.1923-0.1938	0.1935	0.1979-0.1989	0.1981
	1.4	0.3856-0.3872	0.3866	0.3948-0.3958	0.3955
	1.6	0.5812-0.5829	0.5825	0.5940-0.5950	0.5941
1	0.6	0.006136-0.006573	0.006579	0.006307-0.006729	0.006096
	0.8	0.02526-0.02622	0.02581	0.02617-0.02711	0.02520
	1	0.08225-0.08395	0.08250	0.08480-0.08640	0.08593
	1.2	0.1965-0.1989	0.1977	0.2097-0.2115	0.2132
	1.4	0.3560-0.3589	0.3574	0.3946-0.3957	0.3954
	1.6	0.5341-0.5385	0.5364	0.5930-0.5941	0.5932
2	0.6	0.01444-0.01513	0.01529	0.01528-0.01594	0.01365
	0.8	0.04522-0.04655	0.04613	0.04596-0.04719	0.04659
	1	0.1046-0.1067	0.1077	0.1149-0.1168	0.1171
	1.2	0.2054-0.2083	0.2065	0.2319-0.2341	0.2345
	1.4	0.3351-0.3386	0.3382	0.3968-0.3987	0.3991
	1.6	0.4904-0.4944	0.4919	0.5927-0.5938	0.5935

We present the results in Table 1. The option value for both the Bermudan options as well as the European options appears to be accurate. Since the COS method has a very quick convergence, already for $N=64$ the error becomes stable. For at-the-money strikes we have $\log_{10}|\textnormal{error}|\approx 3.5$ . The use of the second-order approximation of the characteristic function is justified by the fact that the option value (and thus the error) stabilizes starting from the second-order approximation. Furthermore, it is noteworthy that the 0th-order approximation is already very accurate.

The computer used in the experiments has an Intel Core i7 CPU with a 2.2 GHz processor. The CPU time of the calculations depends on the number of exercise dates. Assuming we use the second-order approximation of the characteristic function, if we have $M$ exercise dates the CPU time will be $5\cdot M$ ms.

Remark 5.15.

The method can be extended to include time-dependent coefficients. The accuracy and speed of the method will be of the same order as for time-independent coefficients.

Remark 5.16.

The Greeks can be calculated at almost no additional cost using the formulas presented in 3.12. Numerically, the order of convergence is algebraic and is the same for both the exact characteristic function as for the 2nd-order approximation.

5.2 Tests under the CEV-Variance-Gamma dynamics

Consider the jump process to be a Variance-Gamma process. The VG process, is obtained by replacing the time in a Brownian motion with drift $\theta$ and standard deviation $\varrho$ , by a Gamma process with variance $\kappa$ and unitary mean. The model parameters $\varrho$ and $\kappa$ allow to control the skewness and the kurtosis of the distribution of stock price returns. The VG density is characterized by a fat tail and is thus used as a model in situations where small and large asset values are more probable than would be the case for the lognormal distribution. The Lévy measure in this case is given by:

\nu(dx)=\frac{e^{-\lambda_{1}x}}{\kappa x}{\mathds{1}}_{\{x>0\}}dx+\frac{e^{\lambda_{2}x}}{\kappa|x|}{\mathds{1}}_{\{x<0\}}dx,

where

\lambda_{1}=\left(\sqrt{\frac{\theta^{2}\kappa^{2}}{4}+\frac{\varrho^{2}\kappa}{2}}+\frac{\theta\kappa}{2}\right)^{-1},\;\;\;\;\;\lambda_{2}=\left(\sqrt{\frac{\theta^{2}\kappa^{2}}{4}+\frac{\varrho^{2}\kappa}{2}}-\frac{\theta\kappa}{2}\right)^{-1}.

Furthermore we have

	$\displaystyle a(x)=\frac{\sigma_{0}^{2}e^{2(\beta-1)x}}{2},$		(5.101)
	$\displaystyle\mu(t,x)=r+\frac{1}{\kappa}\log\left(1-\kappa\theta-\frac{\kappa\varrho^{2}}{2}\right)-a(x),$		(5.102)
	$\displaystyle\psi(\xi)=-a_{0}(\xi^{2}+i\xi)+ir\xi+i\frac{1}{\kappa}\log\left(1-\kappa\theta-\frac{\kappa\varrho^{2}}{2}\right)\xi-\frac{1}{\kappa}\log\left(1-i\kappa\theta\xi+\frac{\xi^{2}\kappa\varrho^{2}}{2}\right).$		(5.103)

We use the following parameters $S_{0}=1$ , $r=5\%$ , $\sigma_{0}=20\%$ , $\beta=0.5$ , $\kappa=1$ , $\theta=-50\%$ , $\varrho=20\%$ . The results for the European and Bermudan option are presented in Table 2.

Table 2: Prices for a European and a Bermudan Put option (10 exercise dates, expiry

T=1

) in the CEV-VG model for the 2nd-order approximation of the characteristic function, and a Monte Carlo method.

	European		Bermudan
K	MC 95% c.i.	Value	MC 95% c.i.	Value
0.6	0.03090-0.03732	0.03546	0.03756-0.03876	0.03749
0.8	0.08046-0.08247	0.08029	0.08290-0.08484	0.08395
1	0.1507-0.1531	0.1511	0.1572-0.1600	0.1594
1.2	0.2501-0.2538	0.2522	0.2634-0.2668	0.2685
1.4	0.3831-0.3876	0.3847	0.4073-0.4108	0.4137
1.6	0.5430-0.5479	0.5436	0.5920-0.5938	0.5937

5.3 CEV-like Lévy process with a state-dependent measure and default

In this section we consider a model similar to the one used in [7]. The model is defined with local volatility, local default and a state-dependent Lévy measure as follows:

	$\displaystyle a(x)=\frac{1}{2}(b_{0}^{2}+\epsilon_{1}b_{1}^{2}\eta(x)),$
	$\displaystyle\gamma(x)=c_{0}+\epsilon_{2}c_{1}\eta(x),$
	$\displaystyle\nu(x,dz)=\epsilon_{3}\nu_{N}(dz)+\epsilon_{4}\eta(x)\nu_{N}(dz),$
	$\displaystyle\eta(x)=e^{\beta x}.$		(5.104)

We will consider Gaussian jumps, meaning that

\displaystyle\nu_{N}(dz)=\lambda\frac{1}{\sqrt{2\pi\delta^{2}}}\exp\left(\frac{-(z-m)^{2}}{2\delta^{2}}\right)dz.

(5.105)

The regular CEV model has several shortcomings: the volatility for instance drops to zero as the underlying approaches infinity; also the model does not allow the underlying to experience jumps. This model tries to overcome these shortcomings, while still retaining CEV-like behaviour through $\eta(x)$ . The local volatility function $\sigma(x)$ behaves asymptotically like the CEV model, $\sigma(x)\sim\sqrt{\epsilon_{1}}b_{1}e^{\beta x/2}$ as $x\rightarrow-\infty$ , reflecting the fact that the volatility tends to increase as the asset price drops (the leverage effect). Jumps of size $dz$ arrive with a state-dependent intensity of $\nu(x,dz)$ . Lastly, a default arrives with intensity $\gamma(x)$ . The default function $\gamma(x)$ behaves asymptotically like $\epsilon_{2}c_{1}e^{\beta x}$ as $x\rightarrow-\infty$ , reflecting the fact that a default is more likely to occur when the price goes down.
In Table 3 the results are presented for a model as defined in (5.104) without default, meaning that $c_{0}=c_{1}=0$ and with a state-dependent jump measure, so $\nu(x,dz)=\eta(x)\nu_{N}(dz)$ . In this case we have

\psi(\xi)=ir\xi-a_{0}(\xi^{2}-i\xi)-\lambda\nu_{0}(e^{m+\delta^{2}/2}-1)i\xi+\lambda\nu_{0}(e^{mi\xi-\delta^{2}\xi^{2}/2}-1),

where $a_{0}=\frac{1}{2}b_{1}^{2}e^{\beta\bar{x}}$ and $\nu_{0}(dz)=e^{\beta\bar{x}}\nu_{N}(dz)$ . The other parameters are chosen as: $b_{1}=0.15$ , $b_{0}=0$ , $\beta=-2$ , $\lambda=20\%$ , $\delta=20\%$ , $m=-0.2$ , $S_{0}=1$ , $r=5\%$ , $\epsilon_{1}=1$ , $\epsilon_{3}=0$ , $\epsilon_{4}=1$ , the number of exercise dates is 10 and $T=1$ .

Table 3: Prices for a European and a Bermudan Put option (10 exercise dates, expiry

T=1

) in the CEV-like model with state-dependent measure for the 2nd-order approximation characteristic function, and a Monte Carlo method.

	European		Bermudan
K	MC 95% c.i.	Value	MC 95% c.i.	Value
0.8	0.01025-0.01086	0.009385	0.01068-0.01125	0.01024
1	0.04625-0.04745	0.04817	0.05141-0.05253	0.05488
1.2	0.1563-0.1582	0.1564	0.1942-0.1952	0.1952
1.4	0.3313-0.3334	0.3314	0.3927-0.3934	0.3930
1.6	0.5207-0.5229	0.5218	0.5919-0.5926	0.5920
1.8	0.7103-0.7124	0.7122	0.7906-0.7913	0.7910

From the results for both the European option and the Bermudan option we see that the method performs very accurately, even for deeply in-the-money strikes.

In Table 4 the results are presented for the value of a defaultable Put option. In case of default prior to exercise the Put option payoff is 0, in case of no default the value is $(K-S_{t})^{+}$ , depending on the exercise time. We look at the model as defined in (5.104) with the possibility of default and consider state-independent jumps, meaning that we have $\gamma(x)=\eta(x)$ and $\nu(x,dz)=\nu_{N}(dz)$ . We have

\psi(\xi)=ir\xi-a_{0}(\xi^{2}-i\xi)+\gamma_{0}i\xi-\gamma_{0}-\lambda(e^{m+\delta^{2}/2}-1)i\xi+\lambda(e^{mi\xi-\delta^{2}\xi^{2}/2}-1),

where $a_{0}=\frac{1}{2}b_{1}^{2}e^{\beta\bar{x}}$ and $\gamma_{0}=c_{1}e^{\beta\bar{x}}$ . The other parameters are $b_{0}=0$ , $b_{1}=0.15$ , $\beta=-2$ , $c_{0}=0$ , $c_{1}=0.1$ , $S_{0}=1$ , $r=5\%$ , $\epsilon_{1}=1$ , $\epsilon_{2}=1$ , $\epsilon_{3}=1$ , $\epsilon_{4}=0$ , the number of exercise dates is 10 and $T=1$ .

Table 4: Prices for a European and a Bermudan Put option (10 exercise dates, expiry

T=1

) in the CEV-like model with default for the 2nd-order approximation characteristic function, and a Monte Carlo method.

	European		Bermudan
K	MC 95% c.i.	Value	MC 95% c.i.	Value
0.8	0.002905-0.003175	0.003061	0.005876-0.006245	0.006361
1	0.01845-0.01918	0.01893	0.03419-0.03506	0.03520
1.2	0.08148-0.08296	0.08297	0.1820-0.1827	0.1824
1.4	0.2184-0.2205	0.2173	0.3793-0.3801	0.3792
1.6	0.3867-0.3892	0.3841	0.5752-0.5763	0.5763
1.8	0.5597-0.5638	0.5556	0.7727-0.7739	0.7733

Appendix A Proof of Theorem 4.13

Let $X$ and $\tilde{X}$ be as in (4.94) and (4.95) respectively. We first prove that

E[|X_{t}-\tilde{X}_{t}|^{2}]\leq C\left(\kappa_{2}t^{2}+\kappa_{1}^{2}t^{3}\right),\qquad t\in[0,T],

(1.106)

for some positive constant $C$ that depends only on $T$ , on the Lipschitz constants of the coefficients $\mu$ , $\sigma$ , $\eta$ and on the Lévy measure $\nu$ . Here $\kappa_{1}=-\psi^{\prime}(0)$ and $\kappa_{2}=-\psi^{\prime\prime}(0)$ where $\psi$ in (4.96) is the characteristic exponent of the Lévy process $(\tilde{X}_{t}-x)$ .

Using the Hölder inequality, the Itô isometry (see, for instance, [15]) and the Lipschitz continuity of $\eta$ , $\mu$ and $\sigma$ , the mean squared error is bounded by:

$\displaystyle E\left[\|X_{t}-\tilde{X}_{t}\|^{2}\right]\leq$	$\displaystyle\ 3E\left[\left(\int_{0}^{t}(\mu(X_{s})-\mu(x))ds\right)^{2}\right]+3E\left[\left(\int_{0}^{t}(\sigma(X_{s})-\sigma(x))dW_{s}\right)^{2}\right]$	(1.107)
	$\displaystyle+3E\left[\left(\int_{0}^{t}\int_{\mathbb{R}}(\eta(X_{s-})-\eta(x))zd\tilde{N}(s,dz)\right)^{2}\right]$	(1.108)
$\displaystyle\leq$	$\displaystyle\ C\int_{0}^{t}E\left[\|\tilde{X}_{s}-x\|^{2}\right]ds+C\int_{0}^{t}E\left[\|X_{s}-\tilde{X}_{s}\|^{2}\right]ds,$	(1.109)

where

\displaystyle C=6\left(\left\|\mu^{\prime}\right\|_{\infty}^{2}+\left\|\sigma^{\prime}\right\|_{\infty}^{2}+\left\|\eta^{\prime}\right\|_{\infty}^{2}\int_{\mathbb{R}}z^{2}\nu(dz)\right).

(1.110)

Now we recall the following relationship between the first and second moment and cumulants

\displaystyle E[(\tilde{X}_{s}-x)]=c_{1}(s),\qquad E[(\tilde{X}_{s}-x)^{2}]=c_{2}(s)+c_{1}(s)^{2},

(1.111)

where

c_{n}(s)=\frac{s}{i^{n}}\frac{\partial^{n}\psi(\xi)}{\partial\xi^{n}}\bigg{|}_{\xi=0},

and $\psi(\xi)$ is the characteristic exponent of $(\tilde{X}_{s}-x)$ . Thus we have

\displaystyle E\left[|\tilde{X}_{s}-x|^{2}\right]=\kappa_{2}s+\kappa_{1}^{2}s^{2}.

(1.112)

Plugging (1.112) into (1.109) we get

E[|X_{t}-\tilde{X}_{t}|^{2}]\leq C\left(\frac{\kappa_{2}}{2}t^{2}+\frac{\kappa_{1}^{2}}{3}t^{3}\right)+C\int_{0}^{t}E\left[|X_{s}-\tilde{X}_{s}|^{2}\right]ds,

and therefore estimate (1.106) follows by applying the Gronwall inequality in the form

\displaystyle{\varphi}(t)\leq\alpha(t)+C\int_{0}^{t}{\varphi}(s)ds\ \implies\ {\varphi}(t)\leq\alpha(t)+C\int_{0}^{t}{\alpha}(s)e^{C(t-s)}ds,

(1.113)

that is valid for any $C\geq 0$ and ${\varphi}$ , ${\alpha}$ continuous functions.

From (1.106) and (1.112) we can also deduce that

E\left[\left|X_{t}-x\right|^{2}\right]\leq 2E\left[\big{|}X_{t}-\tilde{X}_{t}\big{|}^{2}\right]+2E\left[\big{|}\tilde{X}_{t}-x\big{|}^{2}\right]\leq C\left(\kappa_{2}t+\kappa_{1}^{2}t^{2}\right),\qquad t\in[0,T].

(1.114)

Moreover, from (1.106) we also get the following error estimate for the expectation of a Lipschitz payoff function $v$ :

\left|E\left[v(X_{t})\right]-E[v(\tilde{X}_{t})]\right|\leq C\sqrt{\kappa_{2}t+\kappa_{1}^{2}t^{2}},\qquad t\in[0,T],

(1.115)

where now $C$ also depends on the Lipschitz constant of $v$ . In particular, taking $v(x)=e^{ix{\xi}}$ , this proves (4.97) for $n=0$ .

Next we prove (4.97) for $n=1$ .

Proceeding as in the proof of Lemma 6.23 in [10] with $u(0,x)=\hat{\Gamma}(0,x;t,{\xi})$ and $\bar{x}=x$ , we find

\displaystyle\hat{\Gamma}(0,x;t,{\xi})-\hat{\Gamma}^{(1)}(0,x;t,{\xi})

\displaystyle=\int_{0}^{t}E\left[(L-L_{0})\hat{G}^{1}(s,X_{s};t,{\xi})+(L-L_{1})\hat{G}^{0}(s,X_{s};t,{\xi})\right]ds,

(1.116)

where the 1st-order approximation is as usual

\displaystyle\hat{\Gamma}^{(1)}(s,X;t,{\xi})=\hat{G}^{0}(s,X;t,{\xi})+\hat{G}^{1}(s,X;t,{\xi}),

(1.117)

with

	$\displaystyle\hat{G}^{0}(s,X;t,{\xi})=e^{iX{\xi}+(t-s)\psi({\xi})},$		(1.118)
	$\displaystyle\hat{G}^{1}(s,X;t,{\xi})=e^{iX{\xi}+(t-s)\psi({\xi})}g_{0}^{(1)}(t-s,{\xi}),$		(1.119)

and $g_{0}^{(1)}$ as in (2.45). Using the Lagrangian remainder of the Taylor expansion, we have

$\displaystyle L-L_{0}$	$\displaystyle=\gamma^{\prime}({\varepsilon}^{\prime})(X-x)(\partial_{X}-1)+a^{\prime}({\varepsilon}^{\prime})(X-x)(\partial_{XX}-\partial_{X})+\eta^{\prime}({\varepsilon}^{\prime})(X-x)\int_{\mathbb{R}}\nu(dz)(e^{z}-1-z)\partial_{X}$	(1.120)
	$\displaystyle+\eta^{\prime}({\varepsilon}^{\prime})(X-x)\int_{\mathbb{R}}\nu(dz)(e^{z\partial_{X}}-1-z\partial_{X}),$	(1.121)
$\displaystyle L-L_{1}$	$\displaystyle=\frac{1}{2}\gamma^{\prime\prime}({\varepsilon}^{\prime\prime})(X-x)^{2}(\partial_{X}-1)+\frac{1}{2}a^{\prime\prime}({\varepsilon}^{\prime\prime})(X-x)^{2}(\partial_{XX}-\partial_{X})$	(1.122)
	$\displaystyle+\frac{1}{2}\eta^{\prime\prime}({\varepsilon}^{\prime\prime})(X-x)^{2}\int_{\mathbb{R}}\nu(dz)(e^{z}-1-z)\partial_{X}+\frac{1}{2}\eta^{\prime\prime}({\varepsilon}^{\prime\prime})(X-x)^{2}\int_{\mathbb{R}}\nu(dz)(e^{z\partial_{X}}-1-z\partial_{X}),$	(1.123)

for some ${\varepsilon}^{\prime},{\varepsilon}^{\prime\prime}\in[x,X]$ . Now, $|\hat{G}^{0}|\leq 1$ because $\hat{G}^{0}$ is the characteristic function of the process $\tilde{X}$ in (4.95); thus, we have

\left|(L-L_{1})\hat{G}^{0}(s,X_{s};t,{\xi})\right|\leq C(1+|{\xi}|^{2})\left|X_{s}-x\right|^{2}.

(1.124)

On the other hand, from (2.45) we have

\left|g_{0}^{(1)}(t-s,{\xi})\right|\leq C(t-s)^{2}\left(1+|{\xi}|^{4}\right),

and therefore we get

\displaystyle\left|(L-L_{0})\hat{G}^{1}(s,X_{s};t,{\xi})\right|\leq C(t-s)^{2}(1+|{\xi}|^{4})\left|X_{s}-x\right|.

(1.125)

So we find

\displaystyle\left|\hat{\Gamma}(0,x;t,{\xi})-\hat{\Gamma}^{(1)}(0,x;t,{\xi})\right|

\displaystyle\leq C(1+|{\xi}|^{4})\int_{0}^{t}\left((t-s)^{2}E\left[\left|X_{s}-x\right|\right]+E\left[\left|X_{s}-x\right|^{2}\right]\right)ds

(1.126)

The thesis then follows from estimate (1.114) and integrating.

Appendix B 2nd-order approximation of the characteristic function

For completeness we present here the formulas of the characteristic function approximation in the general case up to the 2nd-order approximation for a process as in (2.1) with a local-volatility coefficient $a(t,x)$ , a local default intensity $\gamma(t,x)$ and a state-dependent measure $\nu(t,x,dz)$ . We expand the coefficients around $\bar{x}=x$ . This choice of $\bar{x}$ is most common in practice and it simplifies the formulas significantly. We have

$\displaystyle\hat{G}^{(0)}(t,x;T,\xi)$	$\displaystyle=e^{i\xi x+(T-t)\psi(\xi)}$	(2.127)
$\displaystyle\hat{G}^{(1)}(t,x;T,\xi)$	$\displaystyle=\hat{G}^{(0)}(t,x;T,\xi)\bigg{(}\frac{1}{2}i(T-t)^{2}\xi(i+\xi)\alpha_{1}\psi^{\prime}(\xi)+\frac{1}{2}(T-t)^{2}(i+\xi)\gamma_{1}\psi^{\prime}(\xi)$	(2.128)
	$\displaystyle-\frac{1}{2}\int_{\mathbb{R}}\nu_{1}(dz)z(T-t)^{2}\xi\psi^{\prime}(\xi)-\frac{1}{2}\int_{\mathbb{R}}\nu_{1}(dz)(e^{z}-1-z)\xi\psi^{\prime}(\xi)$	(2.129)
	$\displaystyle-\frac{1}{2}\int_{\mathbb{R}}i(e^{iz\xi}-1)(T-t)^{2}\psi^{\prime}(\xi)\bigg{)}$	(2.130)
$\displaystyle\hat{G}^{(2)}(t,x;T,\xi)$	$\displaystyle=\hat{G}^{(0)}(t,x;T,\xi)\big{(}G^{(2)}_{1}(t,x;T,\xi)+G^{(2)}_{2}(t,x;T,\xi)+G^{(2)}_{3}(t,x;T,\xi)$	(2.131)
	$\displaystyle+G^{(2)}_{4}(t,x;T,\xi)+G^{(2)}_{5}(t,x;T,\xi)\big{)},$	(2.132)

where we have defined:

$\displaystyle G^{(2)}_{1}(t,x;T,\xi)$	$\displaystyle=\frac{1}{2}((T-t)^{2}a_{2}\xi(i+\xi)\psi^{\prime\prime}(\xi)-\frac{1}{8}(T-t)^{4}a_{1}^{2}\xi^{2}(i+\xi)^{2}\psi^{\prime}(\xi)^{2}$	(2.133)
	$\displaystyle-\frac{1}{6}(T-t)^{3}\xi(i+\xi)(a_{1}^{2}(i+2\xi)\psi^{\prime}(\xi)-2a_{2}\psi^{\prime}(\xi)^{2}+a_{1}^{2}\xi(i+\xi)\psi^{\prime\prime}(\xi)),$	(2.134)
$\displaystyle G^{(2)}_{2}(t,x;T,\xi)$	$\displaystyle=\frac{1}{8}(T-t)^{2}(i+\xi)^{2}\gamma_{1}^{2}\psi^{\prime}(\xi)^{2}+\frac{1}{2}(T-t)^{2}(1-i\xi)\gamma_{2}\psi^{\prime\prime}(\xi)$	(2.135)
	$\displaystyle+\frac{1}{6}(T-t)^{3}(i+\xi)(\gamma_{1}^{2}\psi^{\prime}(\xi)-2i\gamma_{2}\psi^{\prime}(\xi)^{2}+(i+\xi)\gamma_{1}^{2}\psi^{\prime\prime}(\xi)),$	(2.136)
$\displaystyle G^{(2)}_{3}(t,x;T,\xi)$	$\displaystyle=\frac{1}{6}(T-t)^{3}\xi\psi^{\prime}(\xi)\int_{\mathbb{R}^{2}}z\nu_{1}(dz)+\frac{1}{3}i(T-t)^{3}\xi\psi^{\prime}(\xi)^{2}\int_{\mathbb{R}}z\nu_{1}(dz)$	(2.137)
	$\displaystyle+\frac{1}{8}(T-t)^{4}\xi^{2}\psi^{\prime}(\xi)^{2}\int_{\mathbb{R}^{2}}z\nu_{1}(dz)+\frac{1}{2}i\xi(T-t)^{2}\psi^{\prime\prime}(\xi)\int_{\mathbb{R}}z\nu_{1}(dz)$	(2.138)
	$\displaystyle+\frac{1}{6}(T-t)^{3}\xi^{2}\psi^{\prime\prime}(\xi)\int_{\mathbb{R}^{2}}z\nu_{1}(dz),$	(2.139)
$\displaystyle G^{(2)}_{4}(t,x;T,\xi)$	$\displaystyle=-\frac{1}{6}i(T-t)^{3}\psi^{\prime}(\xi)\int_{\mathbb{R}}(e^{iz\xi}-1)\nu_{1}(dz)\int_{\mathbb{R}}ze^{iz\xi}\nu_{1}(dz)$	(2.140)
	$\displaystyle-\frac{1}{8}(T-t)^{4}\psi^{\prime}(\xi)^{2}\int_{\mathbb{R}^{2}}(e^{iz\xi}-1)\nu_{1}(dz)-\frac{1}{3}(T-t)^{3}\psi^{\prime}(\xi)\int_{\mathbb{R}}(e^{iz\xi}-1)\nu_{2}(dz)$	(2.141)
	$\displaystyle-\frac{1}{6}(T-t)^{3}\psi^{\prime\prime}(\xi)\int_{\mathbb{R}^{2}}(e^{iz\xi}-1)\nu_{1}(dz)-\frac{1}{2}(T-t)^{2}\psi^{\prime\prime}(\xi)\int_{\mathbb{R}}(e^{iz\xi}-1)\nu_{2}(dz),$	(2.142)
$\displaystyle G^{(2)}_{5}(t,x;T,\xi)$	$\displaystyle=\frac{1}{6}(T-t)^{3}\xi\psi^{\prime}(\xi)\int_{\mathbb{R}^{2}}(e^{z}-1-z)\nu_{1}(dz)+\frac{1}{8}(T-t)^{4}\xi^{2}\psi^{\prime}(\xi)^{2}\int_{\mathbb{R}^{2}}(e^{z}-1-z)\nu_{1}(dz)$	(2.143)
	$\displaystyle+\frac{1}{3}i(T-t)^{3}\xi\psi^{\prime}(\xi)\int_{\mathbb{R}}(e^{z}-1-z)\nu_{2}(dz)+\frac{1}{6}(T-t)^{3}\xi^{2}\psi^{\prime\prime}(\xi)\int_{\mathbb{R}^{2}}(e^{z}-1-z)\nu_{1}(dz)$	(2.144)
	$\displaystyle+\frac{1}{2}i(T-t)^{2}\xi\psi^{\prime\prime}(\xi)\int_{\mathbb{R}}(e^{z}-1-z)\nu_{2}(dz).$	(2.145)

Essentially $G^{(2)}_{1}$ corresponds to the Taylor expansion of the local volatility, $G^{(2)}_{2}$ results from the default function, $G^{(2)}_{3}$ , $G^{(2)}_{4}$ and $G^{(2)}_{5}$ are related to the state-dependent measure.

References

[1] V. Bally and A. Kohatsu-Higa, A probabilistic interpretation of the parametrix method, Ann. Appl. Probab., 25 (2015), pp. 3095–3138.
[2] A. Capponi, S. Pagliarani, and T. Vargiolu, Pricing vulnerable claims in a Lévy-driven model, Finance Stoch., 18 (2014), pp. 755–789.
[3] P. Carr and V. Linetsky, A jump to default extended CEV model: an application of Bessel processes, Finance Stoch., 10 (2006), pp. 303–330.
[4] F. Fang and C. W. Oosterlee, A novel pricing method for European options based on Fourier-cosine series expansions, SIAM J. Sci. Comput., 31 (2008/09), pp. 826–848.
[5] , Pricing early-exercise and discrete barrier options by Fourier-cosine series expansions, Numer. Math., 114 (2009), pp. 27–62.
[6] S. Heston, A closed-form solution for options with stochastic volatility with applications to bond and currency options, Rev. Financ. Stud., 6 (1993), pp. 327–343.
[7] A. Jacquier and M. Lorig, The smile of certain Lévy-type models, SIAM J. Financial Math., 4 (2013), pp. 804–830.
[8] V. Linetsky, Pricing equity derivatives subject to bankruptcy, Math. Finance, 16 (2006), pp. 255–282.
[9] R. Lord, F. Fang, F. Bervoets, and C. W. Oosterlee, A fast and accurate FFT-based method for pricing early-exercise options under Lévy processes, SIAM J. Sci. Comput., 30 (2008), pp. 1678–1705.
[10] M. Lorig, S. Pagliarani, and A. Pascucci, Analytical expansions for parabolic equations, SIAM J. Appl. Math., 75 (2015), pp. 468–491.
[11] , A family of density expansions for Lévy-type processes, Ann. Appl. Probab., 25 (2015), pp. 235–267.
[12] D. Madan and E. Seneta, The variance gamma (VG) model for share market returns, Journal of Business, 63 (1990), pp. 511–524.
[13] R. Merton, Option pricing when underlying stock returns are discontinuous, Journal of financial economics, 3 (1976), pp. 125–144.
[14] S. Pagliarani, A. Pascucci, and C. Riga, Adjoint expansions in local Lévy models, SIAM J. Financial Math., 4 (2013), pp. 265–296.
[15] A. Pascucci, PDE and martingale methods in option pricing, vol. 2 of Bocconi & Springer Series, Springer, Milan; Bocconi University Press, Milan, 2011.
[16] L. von Sydow, L. J. Höök, E. Larsson, and et al., BENCHOP—the BENCHmarking project in option pricing, Int. J. Comput. Math., 92 (2015), pp. 2361–2379.
[17] B. Zhang and C. W. Oosterlee, Fourier cosine expansions and put-call relations for Bermudan options, in Numerical methods in finance, vol. 12 of Springer Proc. Math., Springer, Heidelberg, 2012, pp. 323–350.