Deng Ding

*Department of Mathematics, University of Macau, Macao, China* 

#### **1. Introduction**

16 Will-be-set-by-IN-TECH

264 Fourier Transform Applications

[1] Hashimoto, I. & Matsumura, A. (2007). Large time behavior of solutions to an initial boundary value problem on the half line for scalar viscous conservation law, *Methods*

[2] Hashimoto, I. & Ueda, Y. (2011). Anti-derivative method in the half space and application to damped wave equations with non-convex convection, to appear in

[3] Il'in, A. M. & Oleinik, O. A. (1964). Behavior of the solution of the Cauchy problem for certain quasilinear equations for unbounded increase of the time, *Amer. Math. Soc.*

[4] Kawashima, S. & Matsumura, A. (1985). Asymptotic stability of traveling wave solutions of systems for one-dimensional gas motion, *Commun. Math. Phys.*, Vol. 101,

[5] Kawashima, S.; Nishibata, S. & Nishikawa, M. (2003). Asymptotic stability of stationary waves for two-dimensional viscous conservation laws in half plane, *Discrete Contin.*

[6] Kawashima, S.; Nishibata, S. & Nishikawa, M. (2004). *L<sup>p</sup>* energy method for multi-dimensional viscous conservation laws and application to the stability of planar

[7] Kawashima, S.; Nishibata, S. & Nishikawa, M. (2004). Asymptotic stability of stationary waves for multi-dimensional viscous conservation laws in half space, preprint. [8] Liu T.-P.; Matsumura, A. & Nishihara, K. (1998). Behaviors of solutions for the Burgers equation with boundary corresponding to rarefaction waves, *SIAM J. Math. Anal.*, Vol.

[9] Liu, T.-P. & Nishihara, K. (1997). Asymptotic behavior for scalar viscous conservation

[10] Nishikawa, M. (1998). Convergence rate to the traveling wave for viscous conservation

[11] Ueda, Y. (2008). Asymptotic stability of stationary waves for damped wave equations with a nonlinear convection term, *Adv. Math. Sci. Appl.*, Vol. 18, No. 1, 329–343. [12] Ueda, Y.; Nakamura, T. & Kawashima, S. (2008). Stability of planar waves for damped wave equations with nonlinear convection in multi-dimensional half space, *Kinetic and*

[13] Ueda, Y.; Nakamura, T. & Kawashima, S. (2010). Stability of degenerate stationary

[14] Ueda, Y.; Nakamura, T. & Kawashima, S. (2011). Energy method in the partial Fourier space and application to stability problems in the half space, *J. Differential Equations*, Vol.

laws with boundary effect, *J. Differential Equations*, Vol. 133, 296–320.

waves for viscous gases, *Arch. Rational Mech. Anal.*, Vol. 198, 735–762.

waves, *J. Hyperbolic Differential Equations*, Vol. 1, 581–603.

**4. References**

*Appl. Anal.*, Vol. 14, 45–60.

*Transl.*, Vol. 42, 19–23.

*Dyn. Syst., Suppl.*, 469–476.

laws, *Funkcial. Ekvac.*, Vol. 41, 107–132.

*Related Models*, Vol. 1, 49–64.

97–127.

29, 293–308.

250, 1169-1199.

*Kyushu Journal of Mathematics*.

Since Bachelier, a French mathematician, first tried to give a mathematical definition for Brownian motion and used it to model the dynamics of stock process in 1900, financial mathematics has developed a lot. Black & Scholes (1973) and Merton (1973) respectively used the geometric Brownian motion (GBM) to model the underlying asset's price process so that they opened the gate of easy ways to compute option prices, which led one of the major breakthroughs of modern finance. Many good ideas have been proposed to model the stock pricing processes since then. Merton (1976) first introduced the jumps into the asset price processes in his seminar paper. More recently, a lot of exponential Lévy models, including Kou's model, Variance Gamma (VG) model, Inverse Gaussian (IG) model, Normal Inverse Gaussian (NIG) model, and CGMY model, etc., were proposed to add jumps in the financial models so that they can describe the statistical properties of financial time series better [e.g. see Cont & Tankov (2003) and references there in]. Also, serval stochastic volatility modes were presented [e.g. see Heston (1993), Bates (1996), and Duffie *et al* (2000), etc.]. Empirical financial data indicate that these models are usually more consist with financial markets.

Under assuming that the price of an underlying asset follows a GBM, Black & Scholes (1973) showed that the value of a European option satisfies a boundary problem of heat equation (Black-Scholes equation) so that they derived an explicit formula (Black-Scholes formula) for the value. However, this nice analytic tractability in option pricing can not be carried over to the most exponential Lévy models or the stochastic volatility modes for asset returns. Thus, many new approaches, including efficient numerical methods, for option pricing have been proposed. These approaches and methods can be classified into four major groups: The partial integro-differential equation (PIDE) and various numerical methods for such equations; The Monte Carlo methods via the stochastic simulation techniques for underlying asset price processes; Directly numerical integration and various numerical methods via differential integral transforms; Backward stochastic differential equation and its numerical methods. Each of them has its advantages and disadvantages for different financial models and specific applications.

Since Stein & Stein (1991) first used Fourier inversion method to find the distribution of the underlying asset in a stochastic volatility model, the Fourier transform methods have become a very active field of financial mathematics. Heston (1993) applied the characteristic function

<sup>\*</sup>The work was supported by the research grant MYRG136(Y1-L2)-FST11-DD from the Research Committee of University of Macau.

where i <sup>=</sup> √−1 is the imaginary unit. Given the function <sup>F</sup>*g*(*u*), the function *<sup>g</sup>*(*x*) can be

Fourier Transform Methods for Option Pricing 267

Sometime it may be more convenient to consider the *generalized Fourier transform* on the

<sup>−</sup>*ax*|*g*(*x*)|*dx* <sup>&</sup>lt; <sup>∞</sup>

<sup>−</sup>*bx*|*g*(*x*)|*dx* <sup>&</sup>lt; <sup>∞</sup>.

<sup>i</sup>*zx <sup>g</sup>*(*x*)*dx*, *<sup>z</sup>* <sup>∈</sup> **<sup>C</sup>**,

*z* ∈ **C** : *a* < Im{*z*} < *b*

*f*(*y*)*g*(*x* − *y*)*dy*.

*f*(*x*) *u* − *x*

*f*(*x*)*g*(*x*)*dx*.

*dx*.

<sup>2</sup>*<sup>π</sup>* �F *f*(*u*), F*g*(*u*)�, where �·, ·� is the inner product of two

<sup>−</sup>i*ux*F*g*(*u*)*du*, *<sup>x</sup>* <sup>∈</sup> **<sup>R</sup>**. (2)

<sup>−</sup>i*xz*F*g*(*z*)*dz*, *<sup>x</sup>* <sup>∈</sup> **<sup>R</sup>**, (3)

(*x*) is the derivative of *g*(*x*).

. Moreover, within this

recovered by the *Fourier inversion formula*:

complex plane. Let *a*, *b* ∈ **R** with *a* < *b*. Assume

Then, the generalized Fourier transform:

paralleled to the real axis:

P1.Differentiation: <sup>F</sup>

which is defined by

P5.Parseval relation: �*<sup>f</sup>* , *<sup>g</sup>*� <sup>=</sup> <sup>1</sup>

square-integrable function defined by

P2.Modulation: <sup>F</sup>

integral:

exists and is analytic for all *z* in the strip **S** =

*g*� (*x*) 

*eλ<sup>x</sup> g*(*x*) 

and

*<sup>g</sup>*(*x*) = <sup>1</sup>

2*π*

 ∞ −∞ *e*

 ∞ −∞ *e*

<sup>F</sup>*g*(*z*) = <sup>∞</sup>

*<sup>g</sup>*(*x*) = <sup>1</sup>

2*π*

−∞ *e*

 i*w*+∞ i*w*−∞

(*u*) = −i*u*F*g*(*u*), where *g*�

(*<sup>f</sup>* <sup>∗</sup> *<sup>g</sup>*)(*x*) = <sup>∞</sup>

<sup>H</sup> *<sup>f</sup>*(*u*) = <sup>1</sup>

�*f* , *g*� =

(*u*) = F*g*(*u* − i*λ*), where *λ* is a constant. P3.Convolution: F(*f* ∗ *g*)(*u*) = F *f*(*u*)F*g*(*u*), where *f* ∗ *g* is the convolution of *f*(*x*) and *g*(*x*),

−∞

P4.Relation with Hilbert transform: F(sgn · *g*)(*u*) = iH(F*g*)(*u*), where sgn(*x*) is the signum function, and H is the Hilbert transform, which is defined by the Cauchy principal value

> *π P*.*V*. ∞ −∞

> > ∞ −∞

with *a* < *w* < *b*. Here, Im{·} denotes taking the imaginary part of argument. The following properties of Fourier transform are useful in this chapter.

strip the generalized Fourier transform may be inverted by integrating along a straight line

*e*

 ∞ −∞ *e*

approach to obtain an analytic representation for the valuation of European options in the Fourier domain. Duffie *et al* (2000) offered a comprehensive survey that Fourier transform methods are applicable to a wide range of stochastic processes, the class of exponential affine diffusions [e.g. see Kwok *et al* (2010) and Schmelzle (2010)].

Carr & Madan (1999) pioneered the use of the fast Fourier transform (FFT) technique by mapping the Fourier transform directly to option prices via the characteristic function. Since then, many efficient numerical methods by using FFT techniques have been proposed, and many authors have discussed these methods in rigorous detail. Lee (2004) extended their method significantly and proved an error analysis for these FFT methods. Lewis (2001) generalized this approach to several general payoff functions via the convolution of generalized Fourier transforms.

Recently, some new ideas to improve the Fourier transform methods have been raised. Fang & Oosterlee (2008) proposed the COS method which is based on the Fourier and Fourier-cosine expansion. Feng & Linetsky (2008) presented a new method which involves the relation between Fourier transform and the Hilbert transform, and the Sinc expansion in Hardy spaces. Lord *et al* (2008) extended the FFT-based methods to the CONV method, which is based on the a quadrature technique and relies heavily on Fourier transform.

In this chapter, we demonstrate the application of Fourier transform as a very effective tool in option pricing theory. Together with the fast Fourier transform technique and other important properties of Fourier transform, we survey serval different methods for pricing European options and some path dependent options under different financial models.

This chapter is organized as follows. In the next section, the mathematical formulations that will be used in this chapter are reviewed. They include a brief discussion on Fourier transform with its important properties; an introduction of discrete Fourier transform and the FFT idea; a definition of characteristic functions; and a brief review on the exponential Lévy models and stochastic volatility models in financial mathematics. In Section 3, several Fourier transform methods for pricing European options are introduced. They include the Black-Scholes type formulas, the FFT methods for signal underlying asset and for multi underlying assets, and the Fourier expansion methods. In Section 4, three new Fourier transform methods for pricing path dependent options are considered. They are the CONV method for pricing Bermudan barrier option, the COS method for pricing Bermudan barrier option, and the fast Hilbert transform method for pricing barrier option.

## **2. Fourier transforms and characteristic functions**

#### **2.1 Fourier transforms**

Let *g*(*x*) be a piecewise continuous real function over **R** which satisfies the integrability condition:

$$\int\_{-\infty}^{\infty} |g(x)| dx < \infty.$$

The *Fourier transform* of *g*(*x*) is defined by

$$\mathcal{F}\mathcal{g}(u) = \int\_{-\infty}^{\infty} e^{iux} \mathcal{g}(\mathbf{x}) d\mathbf{x}, \quad u \in \mathbb{R}, \tag{1}$$

where i <sup>=</sup> √−1 is the imaginary unit. Given the function <sup>F</sup>*g*(*u*), the function *<sup>g</sup>*(*x*) can be recovered by the *Fourier inversion formula*:

$$g(\mathbf{x}) = \frac{1}{2\pi} \int\_{-\infty}^{\infty} e^{-i\mathbf{u}\cdot\mathbf{x}} \mathcal{F}\mathcal{g}(u) du, \quad \mathbf{x} \in \mathbb{R}. \tag{2}$$

Sometime it may be more convenient to consider the *generalized Fourier transform* on the complex plane. Let *a*, *b* ∈ **R** with *a* < *b*. Assume

$$\int\_{-\infty}^{\infty} e^{-ax} |g(x)| dx < \infty$$

and

2 Will-be-set-by-IN-TECH

approach to obtain an analytic representation for the valuation of European options in the Fourier domain. Duffie *et al* (2000) offered a comprehensive survey that Fourier transform methods are applicable to a wide range of stochastic processes, the class of exponential affine

Carr & Madan (1999) pioneered the use of the fast Fourier transform (FFT) technique by mapping the Fourier transform directly to option prices via the characteristic function. Since then, many efficient numerical methods by using FFT techniques have been proposed, and many authors have discussed these methods in rigorous detail. Lee (2004) extended their method significantly and proved an error analysis for these FFT methods. Lewis (2001) generalized this approach to several general payoff functions via the convolution of

Recently, some new ideas to improve the Fourier transform methods have been raised. Fang & Oosterlee (2008) proposed the COS method which is based on the Fourier and Fourier-cosine expansion. Feng & Linetsky (2008) presented a new method which involves the relation between Fourier transform and the Hilbert transform, and the Sinc expansion in Hardy spaces. Lord *et al* (2008) extended the FFT-based methods to the CONV method, which is based on

In this chapter, we demonstrate the application of Fourier transform as a very effective tool in option pricing theory. Together with the fast Fourier transform technique and other important properties of Fourier transform, we survey serval different methods for pricing European

This chapter is organized as follows. In the next section, the mathematical formulations that will be used in this chapter are reviewed. They include a brief discussion on Fourier transform with its important properties; an introduction of discrete Fourier transform and the FFT idea; a definition of characteristic functions; and a brief review on the exponential Lévy models and stochastic volatility models in financial mathematics. In Section 3, several Fourier transform methods for pricing European options are introduced. They include the Black-Scholes type formulas, the FFT methods for signal underlying asset and for multi underlying assets, and the Fourier expansion methods. In Section 4, three new Fourier transform methods for pricing path dependent options are considered. They are the CONV method for pricing Bermudan barrier option, the COS method for pricing Bermudan barrier option, and the fast Hilbert

Let *g*(*x*) be a piecewise continuous real function over **R** which satisfies the integrability


<sup>i</sup>*ux <sup>g</sup>*(*x*)*dx*, *<sup>u</sup>* <sup>∈</sup> **<sup>R</sup>**, (1)

 ∞ −∞

> ∞ −∞ *e*

F*g*(*u*) =

diffusions [e.g. see Kwok *et al* (2010) and Schmelzle (2010)].

the a quadrature technique and relies heavily on Fourier transform.

options and some path dependent options under different financial models.

generalized Fourier transforms.

transform method for pricing barrier option.

The *Fourier transform* of *g*(*x*) is defined by

**2.1 Fourier transforms**

condition:

**2. Fourier transforms and characteristic functions**

$$\int\_{-\infty}^{\infty} e^{-bx} |g(x)| dx < \infty.$$

Then, the generalized Fourier transform:

$$\mathcal{F}\mathfrak{g}(z) = \int\_{-\infty}^{\infty} e^{izx} \mathfrak{g}(x) dx, \quad z \in \mathbb{C}\_{\prime}$$

exists and is analytic for all *z* in the strip **S** = *z* ∈ **C** : *a* < Im{*z*} < *b* . Moreover, within this strip the generalized Fourier transform may be inverted by integrating along a straight line paralleled to the real axis:

$$g(\mathbf{x}) = \frac{1}{2\pi} \int\_{\text{i}w - \infty}^{\text{i}w + \infty} e^{-\text{i}x\overline{z}} \mathcal{F}g(z)dz, \quad \mathbf{x} \in \mathbb{R}, \tag{3}$$

with *a* < *w* < *b*. Here, Im{·} denotes taking the imaginary part of argument.

The following properties of Fourier transform are useful in this chapter.


$$(f \ast g)(x) = \int\_{-\infty}^{\infty} f(y)g(x - y) dy.$$

P4.Relation with Hilbert transform: F(sgn · *g*)(*u*) = iH(F*g*)(*u*), where sgn(*x*) is the signum function, and H is the Hilbert transform, which is defined by the Cauchy principal value integral:

$$\mathcal{H}f(\boldsymbol{\mu}) = \frac{1}{\pi} \boldsymbol{P}.\\ \boldsymbol{V}. \int\_{-\infty}^{\infty} \frac{f(\boldsymbol{x})}{\boldsymbol{\mu} - \boldsymbol{x}} d\boldsymbol{x}.$$

P5.Parseval relation: �*<sup>f</sup>* , *<sup>g</sup>*� <sup>=</sup> <sup>1</sup> <sup>2</sup>*<sup>π</sup>* �F *f*(*u*), F*g*(*u*)�, where �·, ·� is the inner product of two square-integrable function defined by

$$
\langle f, g \rangle = \int\_{-\infty}^{\infty} f(x)g(x) dx.
$$

i.e. it is the Fourier transform of the density *f*(*X*). Here **E**[·] denotes the mathematical expectation. The characteristic functions are uniformly continuous and non-negative definite.

Fourier Transform Methods for Option Pricing 269

P8.*X* and *Y* have the same distribution function if and only if they have the same characteristic

An adapted stochastic process *Xt* with *X*<sup>0</sup> = 0 is called a *Lévy process* if it has independent and stationary increments, and it is continuous in probability. Moreover, every Lévy process has a right continuous with left limits (càdlàd) version which is also a Lévy process. We always

Let *Xt* be a Lévy process, and let <sup>B</sup><sup>0</sup> be the set family of all Borel sets *<sup>U</sup>* <sup>⊂</sup> **<sup>R</sup>** whose closure *<sup>U</sup>*¯

clear that *ν*(*dx*) is a *σ*-finite measure on B0, which is called the *Lévy measure* of *Xt*. Then, the

*tψ*(*u*)

Here the triple (*α*, *σ*2, *ν*(*dx*)) is called the Lévy triple of *Xt*. The proof of this Lévy-Khintchine decomposition is based on the famous *Lévy-Itô decomposition*, and the detail can be found in Cont & Tankov (2003). Also a Lévy process is a strong Markov process, and a semimartingale.

*<sup>e</sup><sup>x</sup>* <sup>−</sup> <sup>1</sup> <sup>−</sup> *<sup>x</sup>*1{|*x*|≤1}(*x*)

To ensure positivity as well as the independence and stationarity of log-returns, in financial mathematics, the price process of the underlying asset is usually modeled as the exponential

There are many models in the financing modeling literature simply correspond to different

, *t* > 0, *U* ∈ B0,

[0, 1], *<sup>U</sup>* for any *<sup>U</sup>* ∈ B0. It is

*ν*(*dx*). (6)

, *u* ∈ **R**, (5)

{|*x*|≥1} *<sup>e</sup>xν*(*dx*) <sup>&</sup>lt; <sup>∞</sup>, the exponential *<sup>e</sup>Xt* is a martingale if

*ν*(*dx*) = 0.

*St* <sup>=</sup> *<sup>S</sup>*0*eXt* , *<sup>t</sup>* <sup>≥</sup> 0. (7)

{0}, *<sup>U</sup>* <sup>≡</sup> 0. Then *<sup>μ</sup>*(*dt*, *dx*) is a *<sup>σ</sup>*-finite random measure on <sup>B</sup>(**R**+) × B0,

 *μ* 

<sup>i</sup>*ux* <sup>−</sup> <sup>1</sup> <sup>−</sup> <sup>i</sup>*ux*1{|*x*|≤1}(*x*)

1*U* Δ*Xs* 

Also they possess the following properties:

function.

**2.4 Exponential Lévy models**

work with this càdlàg version.

*μ* 

where the characteristic exponent function:

Furthermore, under the condition:

*<sup>ψ</sup>*(*u*) = <sup>i</sup>*α<sup>u</sup>* <sup>−</sup> <sup>1</sup>

*α* + 1 2 *σ*<sup>2</sup> + **R** 

choices for *σ*<sup>2</sup> and *ν*(*dx*) to the Lévy process *Xt*.

2 *σ*2*u*<sup>2</sup> + **R** *e*

(0, *t*], *U*

which is called the *jump measure* of *Xt*. Set *ν*(*U*) = **E**

= ∑ 0<*s*≤*t*

characteristic function of Lévy process *Xt* has the *Lévy-Khintchine representation*:

*ϕt*(*u*) = *e*

does not contain 0. Set

and denote *μ*

and only if

of a Lévy process *Xt*:

P6. If *X* and *Y* are independent, then *ϕX*+*Y*(*u*) = *ϕX*(*u*)*ϕY*(*u*). P7. If *<sup>a</sup>*, *<sup>b</sup>* <sup>∈</sup> **<sup>R</sup>** and *<sup>Y</sup>* <sup>=</sup> *aX* <sup>+</sup> *<sup>b</sup>*, then *<sup>ϕ</sup>Y*(*u*) = *<sup>e</sup>*i*buϕX*(*au*).

#### **2.2 Fast Fourier transforms**

Let *<sup>N</sup>* <sup>=</sup> <sup>2</sup>*<sup>L</sup>* be a power of 2, and let *<sup>x</sup>* = (*x*0, *<sup>x</sup>*1,..., *xN*−1)*<sup>T</sup>* be a given *<sup>N</sup>*-dimensional vector, where *v<sup>T</sup>* presents the transpose of a vector *v*. The *discrete Fourier transform* of *x* is another *<sup>N</sup>*-dimensional vector <sup>D</sup>*<sup>x</sup>* = (*y*0, *<sup>y</sup>*1,..., *yN*−1)*<sup>T</sup>* is defined by

$$y\_p = \sum\_{j=0}^{N-1} e^{j\frac{2\pi}{N}jp} x\_{j\prime} \quad p = 0, 1, \ldots, N-1. \tag{4}$$

Denote *AN* an *<sup>N</sup>* <sup>×</sup> *<sup>N</sup>* matrix whose (*p*, *<sup>j</sup>*)th entry is given by

$$a\_{p\,j} = e^{\mathbf{i}\frac{2\pi}{N}jp}, \quad j, p = 0, 1, \dots, N-1.$$

Then, the *discrete Fourier transform* of *x* is given by

$$\mathcal{D}\mathfrak{x} = A\_N \mathfrak{x}.$$

It is clear that the computation to find <sup>D</sup>*<sup>x</sup>* requires *<sup>N</sup>*<sup>2</sup> steps. However, the *fast Fourier technique* (FFT) would require only <sup>1</sup> <sup>2</sup> *NL* <sup>=</sup> <sup>1</sup> <sup>2</sup> *N* log2 *N* steps. The idea behind the FFT algorithm is to take advantage of the periodicity property of *N* roots of unity. For the given *N*-dimensional vector *<sup>x</sup>*, denote *<sup>x</sup>*� = (*x*0, *<sup>x</sup>*2,..., *xN*−2)*<sup>T</sup>* and *<sup>x</sup>*�� = (*x*1, *<sup>x</sup>*3,..., *xN*−1)*<sup>T</sup>* two *<sup>N</sup>*/2-dimensional vectors. Then, we find

*y*� = *AMx*� and *y*�� = *AMx*��, where *M* = *N*/2, *y*� = (*y*� <sup>0</sup>, *y*� <sup>1</sup>,..., *y*� *<sup>M</sup>*−1)*<sup>T</sup>* and *<sup>y</sup>*�� = (*y*�� <sup>0</sup> , *y*�� <sup>1</sup> ,..., *y*�� *<sup>M</sup>*−1)*T*, and the matrix *AM* = [*aj k*]*M*×*<sup>M</sup>* given by

$$a\_{p\,j} = e^{\mathrm{i}\,\frac{2\pi}{N}jp}, \quad j, p = 0, 1, \ldots, M-1.$$

It is easy to verify that the vector D*x*, which is defined in (4), now is given by

$$\begin{aligned} y\_p &= y'\_p + e^{\mathrm{i}\frac{2\pi}{N}p} y\_{p'}^{\prime\prime} \quad p = 0, 1, \ldots, M-1, \\ y\_{M+p} &= y'\_p - e^{\mathrm{i}\frac{2\pi}{N}p} y\_{p'}^{\prime\prime} \quad p = 0, 1, \ldots, M-1. \end{aligned}$$

Instead of performing the matrix-vector multiplication *ANx*, we now only need to perform two matrix-vector multiplications *AMx*� and *AMx*�� so that the number of operations is reduced from *N*<sup>2</sup> to 2(*N*/2)<sup>2</sup> = *N*2/2. The same procedure of reducing the length of the vector by half can be applied repeatedly. Using this FFT algorithm, the total number of operations is reduced from *<sup>O</sup>*(*N*2) to *<sup>O</sup>*(*<sup>N</sup>* log2 *<sup>N</sup>*). The similar FFT algorithm also can be used to calculate the *discrete Fourier inversion transform*:

$$x\_p = \sum\_{j=0}^{N-1} e^{-\mathrm{i}\frac{2\pi}{N}jp} y\_{j\prime} \quad p = 0, 1, \ldots, N-1.$$

#### **2.3 Characteristic functions**

Let *X* be a random variable having the density function *f*(*x*). Then, the *characteristic function* of *X* is defined by

$$\varphi\_X(\mu) = \mathbb{E}\left[e^{\mathrm{i}\mu X}\right] = \int\_{-\infty}^{\infty} e^{\mathrm{i}\mu x} f(\mathbf{x}) d\mathbf{x}, \quad \mu \in \mathbb{R}\_{\prime}$$

4 Will-be-set-by-IN-TECH

Let *<sup>N</sup>* <sup>=</sup> <sup>2</sup>*<sup>L</sup>* be a power of 2, and let *<sup>x</sup>* = (*x*0, *<sup>x</sup>*1,..., *xN*−1)*<sup>T</sup>* be a given *<sup>N</sup>*-dimensional vector, where *v<sup>T</sup>* presents the transpose of a vector *v*. The *discrete Fourier transform* of *x* is another

*<sup>N</sup> jp*, *<sup>j</sup>*, *<sup>p</sup>* <sup>=</sup> 0, 1, . . . , *<sup>N</sup>*−1.

D*x* = *ANx*. It is clear that the computation to find <sup>D</sup>*<sup>x</sup>* requires *<sup>N</sup>*<sup>2</sup> steps. However, the *fast Fourier technique*

take advantage of the periodicity property of *N* roots of unity. For the given *N*-dimensional vector *<sup>x</sup>*, denote *<sup>x</sup>*� = (*x*0, *<sup>x</sup>*2,..., *xN*−2)*<sup>T</sup>* and *<sup>x</sup>*�� = (*x*1, *<sup>x</sup>*3,..., *xN*−1)*<sup>T</sup>* two *<sup>N</sup>*/2-dimensional

*y*� = *AMx*� and *y*�� = *AMx*��,

Instead of performing the matrix-vector multiplication *ANx*, we now only need to perform two matrix-vector multiplications *AMx*� and *AMx*�� so that the number of operations is reduced from *N*<sup>2</sup> to 2(*N*/2)<sup>2</sup> = *N*2/2. The same procedure of reducing the length of the vector by half can be applied repeatedly. Using this FFT algorithm, the total number of operations is reduced from *<sup>O</sup>*(*N*2) to *<sup>O</sup>*(*<sup>N</sup>* log2 *<sup>N</sup>*). The similar FFT algorithm also can be

Let *X* be a random variable having the density function *f*(*x*). Then, the *characteristic function*

*<sup>M</sup>*−1)*<sup>T</sup>* and *<sup>y</sup>*�� = (*y*��

*<sup>N</sup> jp*, *<sup>j</sup>*, *<sup>p</sup>* <sup>=</sup> 0, 1, . . . , *<sup>M</sup>*−1.

*<sup>p</sup>* , *p* = 0, 1, . . . , *M*−1,

*<sup>p</sup>* , *p* = 0, 1, . . . , *M*−1.

*<sup>N</sup> jpyj*, *<sup>p</sup>* <sup>=</sup> 0, 1, . . . , *<sup>N</sup>*−1.

<sup>i</sup>*ux <sup>f</sup>*(*x*)*dx*, *<sup>u</sup>* <sup>∈</sup> **<sup>R</sup>**,

*<sup>N</sup> jpxj*, *<sup>p</sup>* <sup>=</sup> 0, 1, . . . , *<sup>N</sup>*−1. (4)

<sup>2</sup> *N* log2 *N* steps. The idea behind the FFT algorithm is to

<sup>0</sup> , *y*��

<sup>1</sup> ,..., *y*��

*<sup>M</sup>*−1)*T*, and the matrix

*<sup>N</sup>*-dimensional vector <sup>D</sup>*<sup>x</sup>* = (*y*0, *<sup>y</sup>*1,..., *yN*−1)*<sup>T</sup>* is defined by

*N*−1 ∑ *j*=0 *e* i <sup>2</sup>*<sup>π</sup>*

*yp* =

Denote *AN* an *<sup>N</sup>* <sup>×</sup> *<sup>N</sup>* matrix whose (*p*, *<sup>j</sup>*)th entry is given by

Then, the *discrete Fourier transform* of *x* is given by

*ap j* = *e*

<sup>2</sup> *NL* <sup>=</sup> <sup>1</sup>

<sup>0</sup>, *y*�

<sup>1</sup>,..., *y*�

i <sup>2</sup>*<sup>π</sup>*

It is easy to verify that the vector D*x*, which is defined in (4), now is given by

*<sup>p</sup>* + *e* i <sup>2</sup>*<sup>π</sup> <sup>N</sup> py*��

*<sup>p</sup>* − *e* i <sup>2</sup>*<sup>π</sup> <sup>N</sup> py*��

*ap j* = *e*

*yp* = *y*�

*yM*+*<sup>p</sup>* = *y*�

used to calculate the *discrete Fourier inversion transform*:

*xp* =

*ϕX*(*u*) = **E**

*N*−1 ∑ *j*=0 *e* <sup>−</sup><sup>i</sup> <sup>2</sup>*<sup>π</sup>*

> *e* i*uX* = ∞ −∞ *e*

i <sup>2</sup>*<sup>π</sup>*

**2.2 Fast Fourier transforms**

(FFT) would require only <sup>1</sup>

where *M* = *N*/2, *y*� = (*y*�

**2.3 Characteristic functions**

of *X* is defined by

*AM* = [*aj k*]*M*×*<sup>M</sup>* given by

vectors. Then, we find

i.e. it is the Fourier transform of the density *f*(*X*). Here **E**[·] denotes the mathematical expectation. The characteristic functions are uniformly continuous and non-negative definite. Also they possess the following properties:

P6. If *X* and *Y* are independent, then *ϕX*+*Y*(*u*) = *ϕX*(*u*)*ϕY*(*u*).

P7. If *<sup>a</sup>*, *<sup>b</sup>* <sup>∈</sup> **<sup>R</sup>** and *<sup>Y</sup>* <sup>=</sup> *aX* <sup>+</sup> *<sup>b</sup>*, then *<sup>ϕ</sup>Y*(*u*) = *<sup>e</sup>*i*buϕX*(*au*).

P8.*X* and *Y* have the same distribution function if and only if they have the same characteristic function.

#### **2.4 Exponential Lévy models**

An adapted stochastic process *Xt* with *X*<sup>0</sup> = 0 is called a *Lévy process* if it has independent and stationary increments, and it is continuous in probability. Moreover, every Lévy process has a right continuous with left limits (càdlàd) version which is also a Lévy process. We always work with this càdlàg version.

Let *Xt* be a Lévy process, and let <sup>B</sup><sup>0</sup> be the set family of all Borel sets *<sup>U</sup>* <sup>⊂</sup> **<sup>R</sup>** whose closure *<sup>U</sup>*¯ does not contain 0. Set

$$\mu\left(\left(0,t\right],\mathcal{U}\right) = \sum\_{00,\ \mathcal{U}\in \mathcal{B}\_0.$$

and denote *μ* {0}, *<sup>U</sup>* <sup>≡</sup> 0. Then *<sup>μ</sup>*(*dt*, *dx*) is a *<sup>σ</sup>*-finite random measure on <sup>B</sup>(**R**+) × B0, which is called the *jump measure* of *Xt*. Set *ν*(*U*) = **E** *μ* [0, 1], *<sup>U</sup>* for any *<sup>U</sup>* ∈ B0. It is clear that *ν*(*dx*) is a *σ*-finite measure on B0, which is called the *Lévy measure* of *Xt*. Then, the characteristic function of Lévy process *Xt* has the *Lévy-Khintchine representation*:

$$\varphi\_t(u) = \mathfrak{e}^{t\psi(u)}, \quad u \in \mathbb{R}, \tag{5}$$

where the characteristic exponent function:

$$\psi(u) = \mathrm{i}\omega u - \frac{1}{2}\sigma^2 u^2 + \int\_{\mathbb{R}} \left( e^{\mathrm{i}ux} - 1 - \mathrm{i}u \mathbf{x} \mathbf{1}\_{\{|x| \le 1\}}(\mathbf{x}) \right) \nu(d\mathbf{x}).\tag{6}$$

Here the triple (*α*, *σ*2, *ν*(*dx*)) is called the Lévy triple of *Xt*. The proof of this Lévy-Khintchine decomposition is based on the famous *Lévy-Itô decomposition*, and the detail can be found in Cont & Tankov (2003). Also a Lévy process is a strong Markov process, and a semimartingale. Furthermore, under the condition: {|*x*|≥1} *<sup>e</sup>xν*(*dx*) <sup>&</sup>lt; <sup>∞</sup>, the exponential *<sup>e</sup>Xt* is a martingale if and only if

$$\alpha + \frac{1}{2}\sigma^2 + \int\_{\mathbb{R}^3} \left(e^{\mathbf{x}} - 1 - \mathfrak{x} \mathbf{1}\_{\{|\mathbf{x}| \le 1\}}(\mathbf{x})\right) \nu(d\mathbf{x}) = 0.$$

To ensure positivity as well as the independence and stationarity of log-returns, in financial mathematics, the price process of the underlying asset is usually modeled as the exponential of a Lévy process *Xt*:

$$S\_l = S\_0 e^{X\_l}, \quad t \ge 0. \tag{7}$$

There are many models in the financing modeling literature simply correspond to different choices for *σ*<sup>2</sup> and *ν*(*dx*) to the Lévy process *Xt*.

• The Finite Moment Log Stable model, which is considered by Carr & Wu (2003): *Xt* is a

Fourier Transform Methods for Option Pricing 271

Although the class of Lévy processes is quite rich, it is sometime insufficient for multi-period financial modeling. Several models combining jumps and stochastic volatility have been

• The Heston's stochastic volatility model [Heston (1993)]: In this model, the price process

motions with correlation *<sup>ρ</sup>* (i.e. �*BS*, *<sup>B</sup>V*�*<sup>t</sup>* <sup>=</sup> *<sup>ρ</sup>t*), *<sup>κ</sup>* is the mean reversion, *<sup>θ</sup>* is the long-run variance level, and *σ* is the volatility-of-volatility parameter. Although *St* is no longer a Lévy process, under the risk neutral probability the characteristic function of the log-price

2*κθ*/*σ*<sup>2</sup> exp

<sup>d</sup>*St* <sup>=</sup> *<sup>μ</sup>St*d*<sup>t</sup>* <sup>+</sup> <sup>√</sup>*Vt St*d*B<sup>S</sup>*

d*Vt* = *κ*(*θ* − *Vt*)d*t* + *σ*

*<sup>σ</sup>*<sup>2</sup> + i*rtu* + i*x*0*u*

*<sup>γ</sup>* sinh *<sup>γ</sup><sup>t</sup>* 2

where *<sup>x</sup>*<sup>0</sup> <sup>=</sup> log *<sup>s</sup>*0, *<sup>r</sup>* is the interest rate, and *<sup>γ</sup>* <sup>=</sup> *σ*2(*u*<sup>2</sup> <sup>+</sup> <sup>i</sup>*u*)+(*<sup>κ</sup>* <sup>−</sup> <sup>i</sup>*ρσu*)2.

<sup>d</sup>*St* <sup>=</sup> *<sup>μ</sup>St*−d*<sup>t</sup>* <sup>+</sup> <sup>√</sup>*Vt St*−d*B<sup>S</sup>*

*ϕt*(*u*) = *ϕ<sup>D</sup>*

*<sup>γ</sup>* sinh *<sup>γ</sup><sup>t</sup>* 2

d*Vt* = *κ*(*θ* − *Vt*)d*t* + *σ*

log(<sup>1</sup> <sup>+</sup> *<sup>k</sup>*) <sup>∼</sup> *<sup>N</sup>*

*<sup>σ</sup>*<sup>2</sup> <sup>+</sup> <sup>i</sup>(*<sup>r</sup>* <sup>−</sup> *<sup>λ</sup>*¯

<sup>2</sup> <sup>+</sup> *<sup>κ</sup>*−i*ρσ<sup>u</sup>*

where the characteristic function of the diffusion part:

*<sup>t</sup>* (*u*) = exp *κθt*(*κ*−i*ρσu*)

 cosh *<sup>γ</sup><sup>t</sup>*

• The Baste's stochastic volatility model [Bates (1996)]: In this model, the price process *St*

compound Poisson process with intensity *λ* and log-normal distribution of jump sizes such

Similarly to the Heston's model, the characteristic function of *Xt* = log *St* is known in

*k*)*tu* + i*x*0*u*

log(1 + ¯

*<sup>t</sup>* (*u*) · *<sup>ϕ</sup><sup>J</sup>*

<sup>i</sup>*ωtu* <sup>−</sup> (i*σu*)*α<sup>t</sup>* sec

*πα* 2 ,

*t* ,

<sup>√</sup>*Vt* <sup>d</sup>*B<sup>V</sup> t* ,

*<sup>t</sup>* and *<sup>B</sup><sup>V</sup>*

<sup>−</sup> (*u*<sup>2</sup> <sup>+</sup> <sup>i</sup>*u*)*v*<sup>0</sup> *γ* coth *<sup>γ</sup><sup>t</sup>*

> <sup>−</sup> (*u*<sup>2</sup> <sup>+</sup> <sup>i</sup>*u*)*v*<sup>0</sup> *γ* coth *<sup>γ</sup><sup>t</sup>*

<sup>2</sup> + *κ* − i*ρσu*

 ,

*<sup>t</sup>* + *St*−d*Zt*,

<sup>√</sup>*Vt* <sup>d</sup>*B<sup>V</sup> t* ,

*<sup>t</sup>* are two standard Brownian motions with correlation *ρ*, and *Zt* is a

*<sup>k</sup>*) <sup>−</sup> <sup>1</sup> 2 *δ*2, *δ*2 .

*<sup>t</sup>*(*u*),

2*κθ*/*σ*<sup>2</sup> exp

*<sup>t</sup>* are two standard Brownian

<sup>2</sup> + *κ* − i*ρσu*

 ,

finite moment log stable process whose characteristic function is given by

*<sup>ϕ</sup>t*(*u*) = exp

*St* is defined by the system of stochastic differential equations (SDEs):

appeared in the literature. We only list two such models here.

with the initial values *S*<sup>0</sup> = *s*<sup>0</sup> and *V*<sup>0</sup> = *v*0, where *B<sup>S</sup>*

<sup>2</sup> <sup>+</sup> *<sup>κ</sup>*−i*ρσ<sup>u</sup>*

and the variance process *Vt* are given by the system of SDEs:

process *Xt* = log *St* is known in closed form:

*<sup>ϕ</sup>t*(*u*) = exp *κθt*(*κ*−i*ρσu*)

 cosh *<sup>γ</sup><sup>t</sup>*

where *B<sup>S</sup>*

closed form:

*ϕD*

*<sup>t</sup>* and *<sup>B</sup><sup>V</sup>*

that if *k* is its jump size then

where *ω*, *α* and *σ* are parameters.

**2.5 Stochastic volatility models**

The cases that *<sup>σ</sup>*<sup>2</sup> � 0 and *<sup>ν</sup>*(*dx*) is a finite measure, i.e. *Xt* is a *Lévy jump diffusion*:

$$X\_t = \mu t + \sigma B\_t + \sum\_{k=1}^{N\_l} \mathfrak{F}\_{k\prime} \quad t \ge 0,$$

where *Bt* is a standard Brownian motion, {*ξk*} is a sequence of independent and identically distributed random variables with common density *f<sup>ξ</sup>* (*x*), and *Nt* is a Poisson process of intensity *λ*, such that *Bt*, {*ξk*} and *Nt* are mutually independent. In these models, the Lévy triple is given by

$$\mathfrak{a} = \mathfrak{\mu} + \lambda \int\_{\{|\mathbf{x}| < 1\}} \mathfrak{x} f\_{\tilde{\xi}}(\mathbf{x}) d\mathbf{x}, \quad \sigma^2 \quad \text{and} \quad \nu(d\mathbf{x}) = \lambda f\_{\tilde{\xi}}(\mathbf{x}) d\mathbf{x}.$$

These models explain the jump part as the market responses to outside news: good news and bad news. The news arrives according to the Poisson process *Nt*, and the price changes in response according to the jump size *ξk*.

• The earliest of these models is due to Merton (1976). In this model, *ξk*s are normally distributed, with mean *μ<sup>J</sup>* and standard deviation *σJ*, so that the characteristic function of *Xt* is given by

$$\varphi\_t(\mu) = \exp\left\{ \mathrm{i}\omega t \mu - \frac{1}{2} \sigma^2 t \mu^2 + \lambda t \left( e^{\mathrm{i}\mu\_f \mu - \frac{1}{2} \sigma\_f^2 \mu^2} - 1 \right) \right\}.$$

• The second one belongs to Kou (2002). In this model, *ξk*s are non-symmetric double exponentially distributed, and the characteristic function of *Xt* is given by

$$\varphi\_t(\mu) = \exp\left\{ \mathrm{i}\omega t \mu - \frac{1}{2} \sigma^2 t \mu^2 + \lambda t \left( \frac{1 - \eta^2}{1 + \eta^2 \mu^2} e^{\mathrm{i}\omega t} - 1 \right) \right\} \prime$$

where *η* and *κ* are parameters.

The cases that *<sup>σ</sup>*<sup>2</sup> <sup>=</sup> 0 and *<sup>ν</sup>*(**<sup>R</sup>** \ {0}) = <sup>∞</sup>, i.e., the models have infinite activity in jumps. We only list some examples here.

• The VG model, which is proposed by Madan *et al* (1998): *Xt* is a variance gamma process with parameters *σ*, *κ* and *θ*, and its characteristic function is given by:

$$\varphi\_t(\mu) = \left(1 - \mathrm{i}\theta\kappa\mu + \frac{1}{2}\sigma^2\kappa\mu^2\right)^{-t/\kappa}.$$

• The NIG model, which is presented by Barndorff-Nielsen (1997): *Xt* is a normal inverse Gaussian process with parameters *α*, *β* and *δ*, and its characteristic function is given by:

$$g\_t(\mu) = \exp\left\{-\delta t \left(\sqrt{\alpha^2 - (\beta + \mathrm{i}\mu)^2} - \sqrt{\alpha^2 - \beta^2}\right)\right\}.$$

• The CGMY model, which is defined by Carr *et al* (2002): *Xt* is a CGMY process whose characteristic function is given by

$$\varphi\_I(\boldsymbol{\mu}) = \exp\left\{ t \mathbf{C} \boldsymbol{\Gamma}(-\mathbf{Y}) \left[ (M - \mathbf{i}\boldsymbol{\mu})^\mathbf{Y} - M^\mathbf{Y} + (\mathbf{G} + \mathbf{i}\boldsymbol{\mu})^\mathbf{Y} - \mathbf{G}^\mathbf{Y} \right] \right\},$$

where the parameters *C*, *G* and *M* are nonnegative, *Y* < 2, and Γ(*x*) is the gamma function.

6 Will-be-set-by-IN-TECH

where *Bt* is a standard Brownian motion, {*ξk*} is a sequence of independent and identically distributed random variables with common density *f<sup>ξ</sup>* (*x*), and *Nt* is a Poisson process of intensity *λ*, such that *Bt*, {*ξk*} and *Nt* are mutually independent. In these models, the Lévy

These models explain the jump part as the market responses to outside news: good news and bad news. The news arrives according to the Poisson process *Nt*, and the price changes in

• The earliest of these models is due to Merton (1976). In this model, *ξk*s are normally distributed, with mean *μ<sup>J</sup>* and standard deviation *σJ*, so that the characteristic function

• The second one belongs to Kou (2002). In this model, *ξk*s are non-symmetric double

The cases that *<sup>σ</sup>*<sup>2</sup> <sup>=</sup> 0 and *<sup>ν</sup>*(**<sup>R</sup>** \ {0}) = <sup>∞</sup>, i.e., the models have infinite activity in jumps. We

• The VG model, which is proposed by Madan *et al* (1998): *Xt* is a variance gamma process

1 − i*θκu* +

• The NIG model, which is presented by Barndorff-Nielsen (1997): *Xt* is a normal inverse Gaussian process with parameters *α*, *β* and *δ*, and its characteristic function is given by:

• The CGMY model, which is defined by Carr *et al* (2002): *Xt* is a CGMY process whose

(*M* − i*u*)

where the parameters *C*, *G* and *M* are nonnegative, *Y* < 2, and Γ(*x*) is the gamma function.

1 2

*<sup>α</sup>*<sup>2</sup> − (*<sup>β</sup>* + <sup>i</sup>*u*)<sup>2</sup> −

*σ*2*tu*<sup>2</sup> + *λt*

*σ*2*tu*<sup>2</sup> + *λt*

*Nt* ∑ *k*=1

*ξk*, *t* ≥ 0,

*x f<sup>ξ</sup>* (*x*)*dx*, *<sup>σ</sup>*<sup>2</sup> and *<sup>ν</sup>*(*dx*) = *<sup>λ</sup> <sup>f</sup><sup>ξ</sup>* (*x*)*dx*.

 *e* <sup>i</sup>*μJu*<sup>−</sup> <sup>1</sup> <sup>2</sup> *<sup>σ</sup>*<sup>2</sup> *J u*2 − 1 .

 <sup>1</sup> <sup>−</sup> *<sup>η</sup>*<sup>2</sup> <sup>1</sup> <sup>+</sup> *<sup>η</sup>*2*u*<sup>2</sup> *<sup>e</sup>*

*<sup>σ</sup>*2*κu*<sup>2</sup><sup>−</sup>*t*/*<sup>κ</sup>*

.

*<sup>Y</sup>* <sup>−</sup> *<sup>M</sup><sup>Y</sup>* + (*<sup>G</sup>* <sup>+</sup> <sup>i</sup>*u*)

*<sup>α</sup>*<sup>2</sup> − *<sup>β</sup>*<sup>2</sup>

 .

*<sup>Y</sup>* <sup>−</sup> *<sup>G</sup>Y*

,

<sup>i</sup>*κ<sup>u</sup>* <sup>−</sup> <sup>1</sup>

 ,

The cases that *<sup>σ</sup>*<sup>2</sup> � 0 and *<sup>ν</sup>*(*dx*) is a finite measure, i.e. *Xt* is a *Lévy jump diffusion*:

*Xt* = *μt* + *σBt* +

triple is given by

of *Xt* is given by

*α* = *μ* + *λ*

response according to the jump size *ξk*.

where *η* and *κ* are parameters.

only list some examples here.

*ϕt*(*u*) = exp

*ϕt*(*u*) = exp

*ϕt*(*u*) = exp

characteristic function is given by

*ϕt*(*u*) = exp

<sup>i</sup>*αtu* <sup>−</sup> <sup>1</sup> 2

exponentially distributed, and the characteristic function of *Xt* is given by

<sup>i</sup>*αtu* <sup>−</sup> <sup>1</sup> 2

with parameters *σ*, *κ* and *θ*, and its characteristic function is given by:

*ϕt*(*u*) =

 − *δt* 

*tC*Γ(−*Y*)

{|*x*|<1}

• The Finite Moment Log Stable model, which is considered by Carr & Wu (2003): *Xt* is a finite moment log stable process whose characteristic function is given by

$$\varphi\_t(u) = \exp\left\{ \mathrm{i}\omega t u - (\mathrm{i}\sigma u)^\alpha t \sec \frac{\pi \alpha}{2} \right\} \prime$$

where *ω*, *α* and *σ* are parameters.

#### **2.5 Stochastic volatility models**

Although the class of Lévy processes is quite rich, it is sometime insufficient for multi-period financial modeling. Several models combining jumps and stochastic volatility have been appeared in the literature. We only list two such models here.

• The Heston's stochastic volatility model [Heston (1993)]: In this model, the price process *St* is defined by the system of stochastic differential equations (SDEs):

$$\begin{cases} \mathbf{dS}\_{l} = \mu \mathbf{S}\_{l} \mathbf{d}t + \sqrt{V\_{l}} \, \mathbf{S}\_{l} \mathbf{d} \mathbf{B}\_{l}^{S} \\ \mathbf{d}V\_{l} = \kappa (\theta - V\_{l}) \mathbf{d}t + \sigma \sqrt{V\_{l}} \, \mathbf{d} \mathbf{B}\_{l}^{V} \end{cases}$$

with the initial values *S*<sup>0</sup> = *s*<sup>0</sup> and *V*<sup>0</sup> = *v*0, where *B<sup>S</sup> <sup>t</sup>* and *<sup>B</sup><sup>V</sup> <sup>t</sup>* are two standard Brownian motions with correlation *<sup>ρ</sup>* (i.e. �*BS*, *<sup>B</sup>V*�*<sup>t</sup>* <sup>=</sup> *<sup>ρ</sup>t*), *<sup>κ</sup>* is the mean reversion, *<sup>θ</sup>* is the long-run variance level, and *σ* is the volatility-of-volatility parameter. Although *St* is no longer a Lévy process, under the risk neutral probability the characteristic function of the log-price process *Xt* = log *St* is known in closed form:

$$\varphi\_t(\boldsymbol{u}) = \frac{\exp\left\{\frac{\kappa\theta t(\kappa - \mathrm{i}\rho\sigma\boldsymbol{u})}{\sigma^2} + \mathrm{i}\mathrm{r}\boldsymbol{u} + \mathrm{i}\boldsymbol{x}\_0\boldsymbol{u}\right\}}{\left(\cosh\frac{\gamma\boldsymbol{t}}{2} + \frac{\kappa - \mathrm{i}\rho\sigma\boldsymbol{u}}{\gamma}\sinh\frac{\gamma\boldsymbol{t}}{2}\right)^{2\kappa\theta/\sigma^2}}\exp\left\{-\frac{(\boldsymbol{u}^2 + \mathrm{i}\boldsymbol{u})\boldsymbol{v}\_0}{\gamma\coth\frac{\gamma\boldsymbol{t}}{2} + \kappa - \mathrm{i}\rho\sigma\boldsymbol{u}}\right\}\boldsymbol{\zeta}$$

where *<sup>x</sup>*<sup>0</sup> <sup>=</sup> log *<sup>s</sup>*0, *<sup>r</sup>* is the interest rate, and *<sup>γ</sup>* <sup>=</sup> *σ*2(*u*<sup>2</sup> <sup>+</sup> <sup>i</sup>*u*)+(*<sup>κ</sup>* <sup>−</sup> <sup>i</sup>*ρσu*)2.

• The Baste's stochastic volatility model [Bates (1996)]: In this model, the price process *St* and the variance process *Vt* are given by the system of SDEs:

$$\begin{cases} \mathbf{dS}\_{l} = \mu \mathbf{S}\_{l-} \mathbf{d}t + \sqrt{V\_{l}} \, \mathbf{S}\_{l-} \mathbf{d} \mathbf{B}\_{l}^{\mathbf{S}} + \mathbf{S}\_{l-} \mathbf{d} \mathbf{Z}\_{l} \, \mathbf{d} \\ \mathbf{d}V\_{l} = \kappa (\theta - V\_{l}) \mathbf{d}t + \sigma \sqrt{V\_{l}} \, \mathbf{d} \mathbf{B}\_{l}^{V} \end{cases}$$

where *B<sup>S</sup> <sup>t</sup>* and *<sup>B</sup><sup>V</sup> <sup>t</sup>* are two standard Brownian motions with correlation *ρ*, and *Zt* is a compound Poisson process with intensity *λ* and log-normal distribution of jump sizes such that if *k* is its jump size then

$$
\log(1+k) \sim N(\log(1+\bar{k}) - \frac{1}{2}\delta^2, \delta^2).
$$

Similarly to the Heston's model, the characteristic function of *Xt* = log *St* is known in closed form:

$$\boldsymbol{\varrho}\_t(\boldsymbol{\mu}) = \boldsymbol{\varrho}\_t^D(\boldsymbol{\mu}) \cdot \boldsymbol{\varrho}\_t^J(\boldsymbol{\mu})\_\prime$$

where the characteristic function of the diffusion part:

$$\log \theta\_t^D(u) = \frac{\exp\left\{\frac{\mathbf{x}\theta t(\mathbf{x} - \mathbf{i}\rho\sigma u)}{\sigma^2} + \mathbf{i}(r - \lambda\overline{k})tu + \mathbf{i}x\_0u\right\}}{\left(\cosh\frac{\gamma t}{2} + \frac{\mathbf{x} - \mathbf{i}\rho\sigma u}{\gamma}\sinh\frac{\gamma t}{2}\right)^{2\kappa\theta/\sigma^2}} \exp\left\{-\frac{(u^2 + \mathbf{i}u)v\_0}{\gamma \coth\frac{\gamma t}{2} + \kappa - \mathbf{i}\rho\sigma u}\right\}\,,$$

Beginning with Heston's work, many authors used the Fourier inversion methods to solve advanced valuation problems [e.g. see Bates (1996), Schmelzle (2010), and references therein]. As mentioned by Carr & Madan (1999), the numerical integrals based on these methods are generally much faster than finite difference solutions to PDEs or PIDEs. Unfortunately, the FFT cannot be used to evaluate these integrals, since the integrands are singular at the required

Fourier Transform Methods for Option Pricing 273

Recently, several authors applied the Property P5 (Parseval relation) to derive the Fourier inversion formulas in option pricing [e.g. see Dufresne *et al* (2009)]. We have well known that, under the risk-neutral probability, the option price *V*(*S*0, *T*) with the terminal payoff *H*(*ST*) is

**1. The method of Carr and Madan**. Carr & Madan (1999) developed a different method designed to use the fast Fourier transform to price options. They introduced a new technique to calculate the Fourier transform of a modified call option price with respect to the logarithmic strike price so that the fast Fourier transform can be applied to calculate the

Consider a European call with the maturity *T* and the strike price *K*, which is written on a stock whose price process is *St* = *ert*+*Xt* , under a risk-neutral probability. Let *fT*(*x*) be the

where *k* = ln *K* is the log strike price. Note that *VT*(*k*) → *S*<sup>0</sup> = 1 as *k* → −∞, and the function *VT*(*k*) is not square-integrable. Thus, we cannot express the Fourier transform in strike in terms of the characteristic function *ϕT*(*u*) of *XT* and then find a range of strikes by Fourier

*<sup>k</sup>*<sup>+</sup> = *e*

*<sup>T</sup>*(*k* : *α*) = *e*

*fT*(*x*)*e*

for some *α* > 0, which is chosen to improve the integrability. Carr & Madan (1999) showed

*<sup>T</sup>*(*<sup>k</sup>* : *<sup>α</sup>*)(*v*) = <sup>∞</sup>

*H*(*x*)*pT*(*x*)*dx* = *e*

<sup>2</sup>*<sup>π</sup>* �F *<sup>H</sup>*(*u*), <sup>F</sup> *pT*(*u*)�.

<sup>−</sup>*rT* <sup>∞</sup> *k*

*<sup>α</sup>kVT*(*k*)

(1+*α*)*xdx* < ∞.

−∞ *e* <sup>i</sup>*vkV*˜

*<sup>T</sup>*(*k* : *α*)*dk*.

(*e<sup>x</sup>* <sup>−</sup> *<sup>e</sup>*

<sup>−</sup>*rT*�*H*(*x*), *pT*(*x*)�,

*<sup>k</sup>*)*fT*(*x*)*dx*, (10)

<sup>−</sup>*rT* <sup>∞</sup> −∞

where *pT*(*x*) is the density function of *ST*. By the Parseval relation we have

*<sup>V</sup>*(*S*0, *<sup>T</sup>*) = *<sup>e</sup>*−*rT*

evaluation point *u* = 0.

*V*(*S*0, *T*) = *e*

<sup>−</sup>*rT***E** *H*(*ST*) = *e*

**3.2 Fast Fourier transform methods - single asset**

density of *XT*. Consider the price of a call option:

<sup>−</sup>*rT***E**

*<sup>e</sup>XT* <sup>−</sup> *<sup>e</sup>*

To obtain a square-integrable function we consider the modified call price: *V*˜

that a sufficient condition for square-integrability of *V*˜(*k*) is given by ∞ −∞

*<sup>ψ</sup>T*(*<sup>v</sup>* : *<sup>α</sup>*) = <sup>F</sup>*V*˜

*VT*(*k*) = *e*

Consider the Fourier transform of *V*˜(*k* : *α*):

given by

integrals.

inversion.

and the characteristic function of the jump part:

$$\varphi\_t^J(u) = \exp\left\{ t\lambda \left( e^{-\frac{1}{2}\delta^2 u^2 + i\left(\log(1+k) - \frac{1}{2}\delta^2\right)u} - 1 \right) \right\}.$$

#### **3. Fourier transform methods for pricing European options**

#### **3.1 Black-Scholes type formulas**

Gurland (1948) used the Convolution property P3 to derive an formula for calculation of the distribution function *F*(*x*) from the characteristic function *ϕ*(*u*) :

$$F(\mathfrak{x}) = \frac{1}{2} + \frac{1}{2\pi} \int\_0^\infty \frac{e^{\mathrm{i}\mathsf{u}\mathsf{x}}\varphi(-\mathfrak{u}) - e^{-\mathrm{i}\mathsf{u}\mathsf{x}}\varphi(\mathfrak{u})}{\mathrm{i}\mathsf{u}} d\mathfrak{u}.$$

Shephard (1991) gave a simple proof of this formula. From this formula by a simple calculation we can get an important formula for Fourier inversion method:

$$\mathbb{P}(X > \mathbf{x}) = 1 - F(\mathbf{x}) = \frac{1}{2} + \frac{1}{\pi} \int\_0^\infty \text{Re}\left(\frac{e^{-\text{i}\mathbf{x}\cdot\mathbf{x}}\varphi(u)}{\text{i}u}\right) du. \tag{8}$$

Here Re{·} denotes taking the real part of argument. Inside the field of finance, this Fourier inversion method was first considered by Stein & Stein (1991) to find the distribution of the underlying asset price in a stochastic volatility model. Heston (1993) applied the formula (8), to obtain a Black-Scholes type formula for pricing a European call option in the stochastic volatility model.

In fact, with help of the risk-neutral valuation, the price of a European call with spot price *S*<sup>0</sup> and strike price *K* is given by

$$\begin{split} V(S\_0, T; K) &= e^{-rT} \mathbb{E} \left[ (S\_T - K)^+ \right] \\ &= \int\_k^\infty e^{-rT + \mathbf{x}} f\_T(\mathbf{x}) d\mathbf{x} - e^{-rT} K \int\_k^\infty f\_T(\mathbf{x}) d\mathbf{x} \ = I\_1 - e^{-rT} K I\_2 \end{split}$$

where *r* is the interest rate, *T* is the maturity, *k* = log *K*, and *fT*(*x*) is the density function of *XT* = log *ST*. It is clear that the second integral is the probability: Π<sup>2</sup> = **P**(*Xt* > *k*). By the Fourier inversion formula (8) we have

$$\Pi\_2 = \frac{1}{2} + \frac{1}{\pi} \int\_0^\infty \text{Re}\left(\frac{e^{-\text{i}\imath k} \varrho\_T(\mu)}{\text{i}\mu}\right) d\mu\,,$$

where *ϕT*(*u*) is the characteristic function of *XT*. By a change of probability measure, we can verify that the second integral is give by *I*<sup>1</sup> = *S*0Π1, where

$$\Pi\_1 = \frac{1}{2} + \frac{1}{\pi} \int\_0^\infty \text{Re}\left(\frac{e^{-\text{i}\imath k}\wp\_T(\imath - \mathbf{i})}{\text{i}\imath\wp\_T(-\mathbf{i})}\right) d\imath\omega\_1$$

Thus, the price of the European call now is given by a Black-Scholes type formula:

$$V(\mathbf{S}\_0, T; K) = \mathbf{S}\_0 \Pi\_1 - \mathbf{K} e^{-r(T-t)} \Pi\_2. \tag{9}$$

8 Will-be-set-by-IN-TECH

Gurland (1948) used the Convolution property P3 to derive an formula for calculation of the

Shephard (1991) gave a simple proof of this formula. From this formula by a simple calculation

2 + 1 *π*

Here Re{·} denotes taking the real part of argument. Inside the field of finance, this Fourier inversion method was first considered by Stein & Stein (1991) to find the distribution of the underlying asset price in a stochastic volatility model. Heston (1993) applied the formula (8), to obtain a Black-Scholes type formula for pricing a European call option in the stochastic

In fact, with help of the risk-neutral valuation, the price of a European call with spot price *S*<sup>0</sup>

where *r* is the interest rate, *T* is the maturity, *k* = log *K*, and *fT*(*x*) is the density function of *XT* = log *ST*. It is clear that the second integral is the probability: Π<sup>2</sup> = **P**(*Xt* > *k*). By the

where *ϕT*(*u*) is the characteristic function of *XT*. By a change of probability measure, we can

 ∞ 0 Re

 ∞ 0 Re

Thus, the price of the European call now is given by a Black-Scholes type formula:

*<sup>V</sup>*(*S*0, *<sup>T</sup>*; *<sup>K</sup>*) = *<sup>S</sup>*0Π<sup>1</sup> <sup>−</sup> *Ke*−*r*(*T*−*t*)

−*rT* ∞ −∞

<sup>−</sup>*rTK*

 ∞ *k*

 *e*−i*ukϕT*(*u*) i*u*

 *<sup>e</sup>*−i*ukϕT*(*<sup>u</sup>* <sup>−</sup> <sup>i</sup>) i*uϕT*(−i)

(*ST* <sup>−</sup> *<sup>K</sup>*)+ <sup>=</sup> *<sup>e</sup>*

<sup>−</sup>*rT*+*<sup>x</sup> fT*(*x*)*dx* <sup>−</sup> *<sup>e</sup>*

*<sup>e</sup>*i*uxϕ*(−*u*) <sup>−</sup> *<sup>e</sup>*−i*uxϕ*(*u*) i*u*

> ∞ 0 Re

*du*.

*<sup>k</sup>*)<sup>+</sup> *fT*(*x*)*dx*

*fT*(*x*)*dx* = *I*<sup>1</sup> − *e*

*du*. (8)

<sup>−</sup>*rTK I*2,

Π2. (9)

 *e*−i*uxϕ*(*u*) i*u*

(*e<sup>x</sup>* <sup>−</sup> *<sup>e</sup>*

 *du*,

> *du*,

 ∞ 0

<sup>=</sup> <sup>1</sup> <sup>−</sup> *<sup>F</sup>*(*x*) = <sup>1</sup>

and the characteristic function of the jump part:

*<sup>t</sup>*(*u*) = exp

**3. Fourier transform methods for pricing European options**

distribution function *F*(*x*) from the characteristic function *ϕ*(*u*) :

2 + 1 2*π*

we can get an important formula for Fourier inversion method:

*<sup>F</sup>*(*x*) = <sup>1</sup>

<sup>−</sup>*rT***E** 

> <sup>Π</sup><sup>2</sup> <sup>=</sup> <sup>1</sup> 2 + 1 *π*

verify that the second integral is give by *I*<sup>1</sup> = *S*0Π1, where

<sup>Π</sup><sup>1</sup> <sup>=</sup> <sup>1</sup> 2 + 1 *π*

 *tλ e* − 1 <sup>2</sup> *<sup>δ</sup>*2*u*<sup>2</sup>+<sup>i</sup> log(1+¯ *<sup>k</sup>*)<sup>−</sup> <sup>1</sup> 2 *δ*2 *<sup>u</sup>* <sup>−</sup> <sup>1</sup> .

*ϕJ*

**3.1 Black-Scholes type formulas**

**P** *X* > *x* 

volatility model.

and strike price *K* is given by

*V*(*S*0, *T*; *K*) = *e*

Fourier inversion formula (8) we have

= ∞ *k e* Beginning with Heston's work, many authors used the Fourier inversion methods to solve advanced valuation problems [e.g. see Bates (1996), Schmelzle (2010), and references therein]. As mentioned by Carr & Madan (1999), the numerical integrals based on these methods are generally much faster than finite difference solutions to PDEs or PIDEs. Unfortunately, the FFT cannot be used to evaluate these integrals, since the integrands are singular at the required evaluation point *u* = 0.

Recently, several authors applied the Property P5 (Parseval relation) to derive the Fourier inversion formulas in option pricing [e.g. see Dufresne *et al* (2009)]. We have well known that, under the risk-neutral probability, the option price *V*(*S*0, *T*) with the terminal payoff *H*(*ST*) is given by

$$V(\mathbf{S}\_0, T) = e^{-rT} \mathbb{E}\left[H(\mathbf{S}\_T)\right] = e^{-rT} \int\_{-\infty}^{\infty} H(\mathbf{x}) p\_T(\mathbf{x}) d\mathbf{x} = e^{-rT} \langle H(\mathbf{x}), p\_T(\mathbf{x})\rangle\_{\mathcal{H}}$$

where *pT*(*x*) is the density function of *ST*. By the Parseval relation we have

$$V(S\_{0\prime}T) = \frac{e^{-rT}}{2\pi} \langle \mathcal{F}H(\mu), \mathcal{F}p\_T(\mu) \rangle.$$

#### **3.2 Fast Fourier transform methods - single asset**

**1. The method of Carr and Madan**. Carr & Madan (1999) developed a different method designed to use the fast Fourier transform to price options. They introduced a new technique to calculate the Fourier transform of a modified call option price with respect to the logarithmic strike price so that the fast Fourier transform can be applied to calculate the integrals.

Consider a European call with the maturity *T* and the strike price *K*, which is written on a stock whose price process is *St* = *ert*+*Xt* , under a risk-neutral probability. Let *fT*(*x*) be the density of *XT*. Consider the price of a call option:

$$V\_T(k) = e^{-rT} \mathbb{E}\left[\left(e^{X\_T} - e^k\right)^+\right] = e^{-rT} \int\_k^\infty (e^x - e^k) f\_T(\mathbf{x}) d\mathbf{x} \tag{10}$$

where *k* = ln *K* is the log strike price. Note that *VT*(*k*) → *S*<sup>0</sup> = 1 as *k* → −∞, and the function *VT*(*k*) is not square-integrable. Thus, we cannot express the Fourier transform in strike in terms of the characteristic function *ϕT*(*u*) of *XT* and then find a range of strikes by Fourier inversion.

To obtain a square-integrable function we consider the modified call price:

$$
\tilde{V}\_T(k:\mathfrak{a}) = e^{\mathfrak{a}k} V\_T(k)
$$

for some *α* > 0, which is chosen to improve the integrability. Carr & Madan (1999) showed that a sufficient condition for square-integrability of *V*˜(*k*) is given by

$$\int\_{-\infty}^{\infty} f\_T(\mathfrak{x}) e^{(1+\alpha)\mathfrak{x}} d\mathfrak{x} < \infty.$$

Consider the Fourier transform of *V*˜(*k* : *α*):

$$\psi\_T(v:\mathfrak{a}) = \mathcal{F}\check{V}\_T(k:\mathfrak{a})(v) = \int\_{-\infty}^{\infty} e^{\mathrm{i}vk} \tilde{V}\_T(k:\mathfrak{a}) dk.$$

**S***<sup>X</sup>* =

that

*z* : *a* < Im{*z*} < *b*

<sup>−</sup>*rT***E** *H e s*+*rT*+*XT*

Fourier transform of *V*(*s*) by the Parseval relation:

with *a* < *w* < *b*, for a range of initial values *s*.

and hence, the option price is given by

call or the cash-secured put, etc.

<sup>F</sup>*h*(*z*) = <sup>∞</sup>

F*V*(*z*) = *ϕT*(−*z*)

*<sup>V</sup>*(*s*) = *<sup>e</sup>ws*+(1−*w*)(*k*−*rT*)

2*π*

<sup>i</sup>*zsV*(*s*)*ds* = *e*

= *e*

F*V*(*z*) = *e*

*<sup>V</sup>*(*s*) = *<sup>e</sup>*−*rT*

−∞ *e* i*xz*

Hence, the generalized Fourier transform of call option price takes the form:

where **S**¯ *<sup>X</sup>* is the complex conjugate set of **S***X*.

*V*(*s*) = *e*

 ∞ −∞ *e*

for all *z* ∈ **S**. Finally we obtain

, and *h*(*x*) = *H*(*ex*+*rT*) is Fourier integrable in some strip **S***<sup>h</sup>* such

*<sup>s</sup>*+*rT*+*x*)*fT*(*x*)*dx*,

*<sup>s</sup>*+*rT*+*x*)*ds*

*<sup>y</sup>*+*rT*)*dy*,

**<sup>S</sup>***<sup>V</sup>* <sup>=</sup> **<sup>S</sup>**¯ *<sup>X</sup>* <sup>∩</sup> **<sup>S</sup>***<sup>h</sup>* �<sup>=</sup> <sup>∅</sup>,

Fourier Transform Methods for Option Pricing 275

Let *s* = log *S*<sup>0</sup> denote the logarithm of current stock value. Then, the option price is given by

where *fT*(*x*) is density function of *XT*. Under the assumption we can compute the generalized

*fT*(*x*)*dx* <sup>∞</sup>

−∞ *e* <sup>i</sup>*zsH*(*e*

> −∞ *e* <sup>i</sup>*zyH*(*e*

*ϕT*(−*z*)F*h*(*z*)*dz*,

*<sup>k</sup>*+*dx* <sup>=</sup> *<sup>e</sup>k*+i*z*(*k*−*rT*)

<sup>i</sup>*z*(i*<sup>z</sup>* <sup>+</sup> <sup>1</sup>) , 1 <sup>&</sup>lt; Im{*z*} <sup>&</sup>lt; <sup>1</sup> <sup>+</sup> *<sup>α</sup>*.

*<sup>ϕ</sup>T*(−*z*)*e*i*u*(*k*−*rT*−*s*) (i*u* − *w*)(1 + i*u* − *w*)

<sup>i</sup>*z*(i*<sup>z</sup>* <sup>+</sup> <sup>1</sup>) .

*du*,

<sup>−</sup>i*zx fT*(*x*)*dx* <sup>∞</sup>

<sup>−</sup>*rT <sup>ϕ</sup>T*(−*z*)F*h*(*z*), *<sup>z</sup>* <sup>∈</sup> **<sup>S</sup>***V*.

<sup>−</sup>*rT* <sup>∞</sup> −∞ *H*(*e*

 = *e*

<sup>−</sup>*rT* <sup>∞</sup> −∞

<sup>−</sup>*rT* <sup>∞</sup> −∞ *e*

Option prices can now be given by the generalized Fourier inversion formula (3):

 i*w*+∞ i*w*−∞

For a European call option, the function *<sup>h</sup>*(*x*)=(*erT*+*<sup>x</sup>* <sup>−</sup> *<sup>e</sup>k*)+, is Fourier integrable in the region Im{*z*} > 1, where its generalized Fourier transform can be computed explicitly:

*<sup>e</sup>x*+*rT* <sup>−</sup> *<sup>e</sup>*

*e*(1+i*z*)(*k*−*rT*)

 ∞ −∞

for some 1 < *w* < 1 + *α*. The integral in this formula can be approximated by using the FFT algorithm as the formula (12). However, the choice of *w* is a delicate issue because choosing big *w* leads to slower decay rates at infinity and bigger truncation errors, and while *w* is close to one the denominator diverges and the discretization error becomes to large (see Chapter 11 in Cont & Tankov (2003)). On the other hand, Lewis (2001) also listed the generalized Fourier transforms F*h*(*z*) and the strip **S***<sup>h</sup>* for various claims, for instance, the put option, the covered

2*π*

By (10) and the definition of characteristic functions we have

$$\psi\_T(v:a) = e^{-rT} \int\_{-\infty}^{\infty} f\_T(\mathbf{x}) d\mathbf{x} \int\_{-\infty}^{\chi} \left( e^{\mathbf{x} + ak} - e^{(1+a)k} \right) e^{i\mathbf{x}k} dk$$

$$= \frac{e^{-rT} \varphi\_T(v - \mathbf{i}(a+1))}{a^2 + a - v^2 + \mathbf{i}(2a+1)v} \tag{11}$$

where *ϕT*(*u*) is the characteristic function of *XT*. Using the Fourier inverse transform we get

$$V\_T(k) = e^{-ak}\tilde{V}(k:\mathfrak{a}) = \frac{e^{-ak}}{2\pi} \int\_{-\infty}^{\infty} e^{-i\upsilon k} \psi\_T(v:\mathfrak{a}) dv.$$

Note that *ψT*(*v* : *α*) is odd in its imaginary part and even in its real part. Since *V*(*k*) is real we have

$$V\_T(k) = \frac{e^{-uk}}{\pi} \int\_0^\infty e^{-i\nu k} \psi\_T(v:a) dv. \tag{12}$$

Now this integral can be evaluated by the numerical approximation using the trapezoidal rule:

$$V\_T(k\_p) \approx \frac{e^{-\alpha k\_p}}{2\pi} \sum\_{j=0}^{N-1} e^{-\text{i}\upsilon\_j k\_p} \psi\_T(\upsilon\_j : \alpha) \Delta \upsilon\_\prime \quad p = 0, 1, \dots, N-1.$$

to apply the FFT algorithm, we set *N* as a power of 2, and define the grid points:

$$\begin{cases} k\_p = -\frac{1}{2}N\Delta k + p\Delta k\_\prime & p = 0, 1, \dots, N-1, \\ v\_j = -\frac{1}{2}N\Delta u + j\Delta v\_\prime & j = 0, 1, \dots, N-1, \end{cases}$$

with the step sizes Δ*k* and Δ*v*, which satisfy the Nyquist relation: Δ*k*Δ*v* = 2*π*/*N*. Then, the numerical approximation can expressed as

$$V\_T(k\_p) \approx \frac{e^{-ak\_p}}{2\pi} \sum\_{j=0}^{N-1} (-1)^{j+p} e^{-i\frac{2\pi}{N}jp} \psi\_T(v\_j : \mathfrak{a}) \Delta v, \quad p = 0, 1, \dots, N-1,$$

which can be efficiently computed by using FFT algorithm.

Many authors have discussed this Carr and Madan's pricing method in rigorous detail. Lee (2004) extended this method significantly and proved an error analysis for this FFT method. On the other hand, numerical experiments show that this FFT method has a quite large error when the strike price *K* is small, and its stability is very dependent on the choice of the damping exponential factor *α*. Ding & U (2010) presented a modified FFT-based method to overcome these disadvantages.

**2. The method of Lewis**. Lewis (2001) introduced an option pricing method which generalizes previous work on Fourier transform methods. The main idea of his method is to express the option price via the convolution of generalized Fourier transforms, and then, to apply the property P5 (Parseval relation) to obtain the generalized Fourier transform of option price.

Consider a European type option whose payoff is *H*(*ST*) at maturity *T*, where *St* = *S*0*ert*+*Xt* is the stock price process under the risk neutral probability. To proceed, we assume that *Xt* is a Lévy process having the analytic characteristic function *ϕT*(*z*), which is regular in the strip 10 Will-be-set-by-IN-TECH

 *x* −∞ 

 ∞ −∞ *e*

*<sup>e</sup>x*+*α<sup>k</sup>* <sup>−</sup> *<sup>e</sup>*

<sup>−</sup>i*vjkpψT*(*vj* : *<sup>α</sup>*)Δ*v*, *<sup>p</sup>* <sup>=</sup> 0, 1, . . . , *<sup>N</sup>* <sup>−</sup> 1.

*<sup>N</sup> jpψT*(*vj* : *<sup>α</sup>*)Δ*v*, *<sup>p</sup>* <sup>=</sup> 0, 1, . . . , *<sup>N</sup>*−1,

<sup>2</sup> *N*Δ*k* + *p*Δ*k*, *p* = 0, 1, . . . , *N*−1,

<sup>2</sup> *N*Δ*u* + *j*Δ*v*, *j* = 0, 1, . . . , *N*−1,

with the step sizes Δ*k* and Δ*v*, which satisfy the Nyquist relation: Δ*k*Δ*v* = 2*π*/*N*. Then, the

Many authors have discussed this Carr and Madan's pricing method in rigorous detail. Lee (2004) extended this method significantly and proved an error analysis for this FFT method. On the other hand, numerical experiments show that this FFT method has a quite large error when the strike price *K* is small, and its stability is very dependent on the choice of the damping exponential factor *α*. Ding & U (2010) presented a modified FFT-based method to

**2. The method of Lewis**. Lewis (2001) introduced an option pricing method which generalizes previous work on Fourier transform methods. The main idea of his method is to express the option price via the convolution of generalized Fourier transforms, and then, to apply the property P5 (Parseval relation) to obtain the generalized Fourier transform of option price. Consider a European type option whose payoff is *H*(*ST*) at maturity *T*, where *St* = *S*0*ert*+*Xt* is the stock price process under the risk neutral probability. To proceed, we assume that *Xt* is a Lévy process having the analytic characteristic function *ϕT*(*z*), which is regular in the strip

<sup>−</sup><sup>i</sup> <sup>2</sup>*<sup>π</sup>*

(1+*α*)*k e* <sup>i</sup>*vkdk*

<sup>−</sup>i*vkψT*(*v* : *α*)*dv*.

, (11)

<sup>−</sup>i*vkψT*(*v* : *α*)*dv*. (12)

*fT*(*x*)*dx*

*<sup>α</sup>*<sup>2</sup> + *<sup>α</sup>* − *<sup>v</sup>*<sup>2</sup> + <sup>i</sup>(2*<sup>α</sup>* + <sup>1</sup>)*<sup>v</sup>*

*v* − i(*α* + 1)

where *ϕT*(*u*) is the characteristic function of *XT*. Using the Fourier inverse transform we get

2*π*

Note that *ψT*(*v* : *α*) is odd in its imaginary part and even in its real part. Since *V*(*k*) is real we

Now this integral can be evaluated by the numerical approximation using the trapezoidal

 ∞ 0 *e*

<sup>−</sup>*<sup>α</sup>kV*˜(*<sup>k</sup>* : *<sup>α</sup>*) = *<sup>e</sup>*−*α<sup>k</sup>*

*π*

to apply the FFT algorithm, we set *N* as a power of 2, and define the grid points:

*VT*(*k*) = *<sup>e</sup>*−*α<sup>k</sup>*

*N*−1 ∑ *j*=0 *e*

By (10) and the definition of characteristic functions we have

−*rT* ∞ −∞

<sup>=</sup> *<sup>e</sup>*−*rT <sup>ϕ</sup><sup>T</sup>*

*ψT*(*v* : *α*) = *e*

*VT*(*k*) = *e*

*VT*(*kp*) <sup>≈</sup> *<sup>e</sup>*−*αkp*

numerical approximation can expressed as

2*π*

*VT*(*kp*) <sup>≈</sup> *<sup>e</sup>*−*αkp*

overcome these disadvantages.

2*π*

*kp* <sup>=</sup> <sup>−</sup><sup>1</sup>

*vj* <sup>=</sup> <sup>−</sup><sup>1</sup>

*N*−1 ∑ *j*=0

which can be efficiently computed by using FFT algorithm.

(−1)*j*<sup>+</sup>*pe*

have

rule:

**S***<sup>X</sup>* = *z* : *a* < Im{*z*} < *b* , and *h*(*x*) = *H*(*ex*+*rT*) is Fourier integrable in some strip **S***<sup>h</sup>* such that

$$\mathbb{S}\_V = \mathsf{\bar{S}}\_X \cap \mathbb{S}\_h \neq \mathcal{O}\_r$$

where **S**¯ *<sup>X</sup>* is the complex conjugate set of **S***X*.

Let *s* = log *S*<sup>0</sup> denote the logarithm of current stock value. Then, the option price is given by

$$V(\mathbf{s}) = e^{-rT} \mathbb{E}\left[ H(e^{\mathbf{s} + rT + X\_T}) \right] = e^{-rT} \int\_{-\infty}^{\infty} H(e^{\mathbf{s} + rT + \mathbf{x}}) f\_T(\mathbf{x}) d\mathbf{x} \lambda$$

where *fT*(*x*) is density function of *XT*. Under the assumption we can compute the generalized Fourier transform of *V*(*s*) by the Parseval relation:

$$\begin{aligned} \int\_{-\infty}^{\infty} e^{\mathrm{i}zs} V(s) ds &= e^{-rT} \int\_{-\infty}^{\infty} f\_T(\mathbf{x}) d\mathbf{x} \int\_{-\infty}^{\infty} e^{\mathrm{i}zs} H(e^{\mathbf{s} + rT + \mathbf{x}}) ds \\ &= e^{-rT} \int\_{-\infty}^{\infty} e^{-\mathrm{i}z\mathbf{x}} f\_T(\mathbf{x}) d\mathbf{x} \int\_{-\infty}^{\infty} e^{\mathrm{i}zy} H(e^{\mathbf{y} + rT}) dy \end{aligned}$$

for all *z* ∈ **S**. Finally we obtain

$$\mathcal{F}V(z) = e^{-rT} \varphi\_T(-z) \mathcal{F}h(z), \quad z \in \mathbb{S}\_V.$$

Option prices can now be given by the generalized Fourier inversion formula (3):

$$V(s) = \frac{e^{-rT}}{2\pi} \int\_{\text{i}w - \infty}^{\text{i}w + \infty} \varphi\_T(-z) \mathcal{F}h(z)dz.$$

with *a* < *w* < *b*, for a range of initial values *s*.

For a European call option, the function *<sup>h</sup>*(*x*)=(*erT*+*<sup>x</sup>* <sup>−</sup> *<sup>e</sup>k*)+, is Fourier integrable in the region Im{*z*} > 1, where its generalized Fourier transform can be computed explicitly:

$$\mathcal{F}h(z) = \int\_{-\infty}^{\infty} e^{i\mathbf{x}z} \left(e^{\mathbf{x}+rT} - e^k\right)^+ d\mathbf{x} = \frac{e^{k+\mathbf{i}z(k-rT)}}{\mathbf{i}z(\mathbf{i}z+1)}.$$

Hence, the generalized Fourier transform of call option price takes the form:

$$\mathcal{F}V(z) = \varphi\_T(-z) \frac{e^{(1+iz)(k-rT)}}{\text{iz}(\text{iz}+1)}, \quad 1 < \text{Im}\{z\} < 1+\mathfrak{a}.$$

and hence, the option price is given by

$$V(s) = \frac{e^{ws + (1 - w)(k - rT)}}{2\pi} \int\_{-\infty}^{\infty} \frac{q\tau\_T(-z)e^{\mathrm{i}\mu(k - rT - s)}}{(\mathrm{i}\mathrm{i} - w)(1 + \mathrm{i}\mathrm{i} - w)} d\mu\_s$$

for some 1 < *w* < 1 + *α*. The integral in this formula can be approximated by using the FFT algorithm as the formula (12). However, the choice of *w* is a delicate issue because choosing big *w* leads to slower decay rates at infinity and bigger truncation errors, and while *w* is close to one the denominator diverges and the discretization error becomes to large (see Chapter 11 in Cont & Tankov (2003)). On the other hand, Lewis (2001) also listed the generalized Fourier transforms F*h*(*z*) and the strip **S***<sup>h</sup>* for various claims, for instance, the put option, the covered call or the cash-secured put, etc.

for *k*<sup>1</sup>

*<sup>p</sup>* � *<sup>k</sup>*<sup>2</sup>

*<sup>q</sup>*, where

*q*)=(−1)*p*+*<sup>q</sup>*

which can be computed via the FFT algorithm.

*N*−1 ∑ *m*=0

the terminal spread option payoff with the strike price *K* = 1, i.e.

<sup>Γ</sup>(*z*) = <sup>∞</sup> 0 *e* −*t t*

> i*�*2+∞ i*�*2−∞

And denote Γ(*z*) the complex gamma function defined by

(2*π*)<sup>2</sup>

*<sup>h</sup>*(*u*1, *<sup>u</sup>*2) = <sup>Γ</sup>

i(*u*1*X*1*T*+*u*2*X*2*<sup>T</sup>* )

*<sup>e</sup>X*1*<sup>T</sup>* <sup>−</sup> *<sup>e</sup>X*2*<sup>T</sup>* <sup>−</sup> <sup>1</sup>

 i*�*2+∞ i*�*2−∞

 i*�*1+∞ i*�*1−∞

 i*�*1+∞ i*�*1−∞

ˆ

function of *X*1*<sup>T</sup>* − *x*<sup>1</sup> and *X*2*<sup>T</sup>* − *x*2. Then, we have **E** *e*

two-dimensional Fourier inversion transform:

= *e* <sup>−</sup>*rT***E**

<sup>=</sup> *<sup>e</sup>*−*rT* (2*π*)<sup>2</sup>

<sup>=</sup> *<sup>e</sup>*−*rT* (2*π*)<sup>2</sup>

<sup>−</sup>*rT***E**

 1 (2*π*)<sup>2</sup>

 i*�*2+∞ i*�*2−∞

 i*�*2+∞ i*�*2−∞

*VT*(*x*1, *x*2) = *e*

*<sup>h</sup>*(*x*1, *<sup>x</sup>*2) = <sup>1</sup>

*N*−1 ∑ *n*=0 *e* <sup>−</sup><sup>i</sup> <sup>2</sup>*<sup>π</sup>*

*<sup>N</sup>* (*mp*+*nq*)

Fourier Transform Methods for Option Pricing 277

There are some other types of terminal payoff functions that admit analytic representation of the Fourier transform of the damped option price [see e.g. Eberlein *et al* (2010), Lee (2004), and references therein]. However, to derive the FFT pricing algorithm for the spread option:

Dempster & Hong (2002) approximated the exercise region of *H*(*S*1*T*, *S*2*T*) by a combination

**2**. Hurd & Zhou (2010) proposed an alternative approach to pricing the European spread option (13) under the three-factor SV model and exponential Lévy models. Let *h*(*x*1, *x*2) be

Then, the method of Hurd and Zhou mainly relies on the following Fourier representation of

 i*�*1+∞ i*�*1−∞

where *�*<sup>1</sup> and *�*<sup>2</sup> are any two given real numbers with *�*<sup>2</sup> > 0 and *�*<sup>1</sup> + *�*<sup>2</sup> < −1, and

 = *e* *e*

i(*u*<sup>1</sup> + *u*2) − 1

The proof of this formula can be found in the appendix of Hurd & Zhou (2010). Let *X*1*<sup>t</sup>* = log *S*1*<sup>t</sup>* and *X*2*<sup>t</sup>* = log *S*2*<sup>t</sup>* with *X*<sup>10</sup> = *x*<sup>1</sup> and *X*<sup>20</sup> = *x*2, and let *ϕT*(*u*1, *u*2) be the characteristic

Now, using the formula (16), the price of the spread option (14) is expressed as an explicit

<sup>+</sup>

 i*�*1+∞ i*�*1−∞

**E** *e*

*e*

*e*

i(*u*1*X*1*T*+*u*2*X*2*<sup>T</sup>* )

<sup>i</sup>(*u*<sup>1</sup> *<sup>x</sup>*1+*u*<sup>2</sup> *<sup>x</sup>*2)*ϕT*(*u*1, *u*2)ˆ

(−1)*m*+*nψT*(*v*<sup>1</sup>

*<sup>H</sup>*(*S*1*T*, *<sup>S</sup>*2*T*)=(*S*1*<sup>T</sup>* <sup>−</sup> *<sup>S</sup>*2*<sup>T</sup>* <sup>−</sup> *<sup>K</sup>*)+, (13)

*<sup>h</sup>*(*x*1, *<sup>x</sup>*2)=(*ex*<sup>1</sup> <sup>−</sup> *<sup>e</sup>x*<sup>2</sup> <sup>−</sup> <sup>1</sup>)+. (14)

*<sup>z</sup>*−1*dt*, Re{*z*} <sup>&</sup>gt; 0.

i(*u*<sup>1</sup> *x*1+*u*<sup>2</sup> *x*2) ˆ

 Γ(−i*u*2)

<sup>Γ</sup>(i*u*<sup>1</sup> <sup>+</sup> <sup>1</sup>) .

<sup>i</sup>(*u*<sup>1</sup> *<sup>x</sup>*1+*u*<sup>2</sup> *<sup>x</sup>*2)*ϕT*(*u*1, *u*2).

i(*u*1*X*1*T*+*u*2*X*2*<sup>T</sup>* ) ˆ

 ˆ *<sup>m</sup>*, *<sup>v</sup>*<sup>2</sup>

*n*)Δ*v*1Δ*v*<sup>2</sup>

*h*(*u*1, *u*2)*du*1*du*2, (15)

*h*(*u*1, *u*2)*du*1*du*<sup>2</sup>

*h*(*u*1, *u*2)*du*1*du*2.

*h*(*u*1, *u*2)*du*1*du*<sup>2</sup>

,

Ψ*T*(*k*<sup>1</sup> *<sup>p</sup>*, *<sup>k</sup>*<sup>2</sup>

of rectangular strips.

the payoff function *h*(*x*1, *x*2):

#### **3.3 Fast Fourier transform methods - multi assets**

**1**. The most direct extension of the Carr and Madan method to the multi assets model is to price the correlation option, whose payoff at the maturity *T* is defined by

$$H(\mathbf{S}\_{1T}, \mathbf{S}\_{2T}) = (\mathbf{S}\_{1T} - \mathbf{K}\_1)^+ (\mathbf{S}\_{2T} - \mathbf{K}\_2)^+ \prime$$

where *K*<sup>1</sup> and *K*<sup>2</sup> are strike prices [see, e.g. Kwok *et al* (2010)]. Define *Xit* = log *Sit* and *ki* = log *Ki*, *i* = 1, 2, and let *fT*(*x*1, *x*2) be the joint density of *X*1*<sup>T</sup>* and *X*2*<sup>T</sup>* under the risk neutral probability. Then, the characteristic function of *X*1*<sup>T</sup>* and *X*2*<sup>T</sup>* is defined by the following two-dimensional Fourier transform:

$$
\varphi\_T(\mu\_1, \mu\_2) = \int\_{-\infty}^{\infty} \int\_{-\infty}^{\infty} e^{\mathbf{i}(\mu\_1 \mathbf{x}\_1 + \mu\_2 \mathbf{x}\_2)} f\_T(\mathbf{x}\_1, \mathbf{x}\_2) d\mathbf{x}\_1 d\mathbf{x}\_2.
$$

Consider the price of the correlation option:

$$V\_T(k\_1, k\_2) = e^{-rT} \int\_{k\_2}^{\infty} \int\_{k\_1}^{\infty} (e^{\mathbf{x}\_1} - e^{k\_1})(e^{\mathbf{x}\_2} - e^{k\_2}) f\_T(\mathbf{x}\_1, \mathbf{x}\_2) d\mathbf{x}\_1 d\mathbf{x}\_2.$$

Following the Carr and Madam method we consider the modified call price:

$$
\tilde{V}\_T(k\_1, k\_2 : \alpha\_1, \alpha\_2) = e^{\alpha\_1 k\_2 + \alpha\_2 k\_2} V\_T(k\_1, k\_2).
$$

for some parameters *α*<sup>1</sup> > 0 and *α*<sup>2</sup> > 0, which are chosen such that this modified price is square integrable for negative value of *k*<sup>1</sup> and *k*2. Then, by a direct integral the Fourier transform of *V*˜(*k*1, *k*2) is given by

$$\begin{split} \psi\_T(v\_1, v\_2) &= \, \_T\tilde{V}\_T(k\_1, k\_2: \mathfrak{a}\_1, \mathfrak{a}\_2) \\ &= \frac{e^{-rT} \varrho\_T \left( v\_1 - \mathrm{i}(\mathfrak{a}\_1 + 1), v\_2 - \mathrm{i}(\mathfrak{a}\_2 + 1) \right)}{(\mathfrak{a}\_1 + \mathrm{i} v\_1)(\mathfrak{a}\_1 + 1 + \mathrm{i} v\_1)(\mathfrak{a}\_2 + \mathrm{i} v\_2)(\mathfrak{a}\_2 + 1 + \mathrm{i} v\_2)}. \end{split}$$

Applying the Fourier inversion on *ψT*(*v*1, *v*2) and using the numerical approximation of two-dimensional Fourier inversion integral we have

$$\begin{split} V\_T(k\_p^1, k\_q^2) &= \frac{\varepsilon^{-a\_1 k\_p^1 - a\_2 k\_q^2}}{(2\pi)^2} \int\_{-\infty}^{\infty} \int\_{-\infty}^{\infty} e^{-\mathbf{i}(v\_1 k\_p^1 + v\_2 k\_q^2)} \psi\_T(v\_1, v\_2) d\upsilon\_1 d\upsilon\_2 \\ &\approx \frac{\varepsilon^{-a\_1 k\_p^1 - a\_2 k\_q^2}}{(2\pi)^2} \sum\_{m=0}^{N-1} \sum\_{n=0}^{N-1} e^{-\mathbf{i}(v\_w^1 k\_p^1 + v\_n^2 k\_q^2)} \psi\_T(v\_m^1, v\_n^2) \Delta v^1 \Delta v^2 \end{split}$$

for *p*, *q* = 0, 1, . . . , *N*−1. To apply the two-dimensional of the FFT algorithm we set

$$v v\_m^1 = \left(m - \frac{N}{2}\right) \Delta v^1, \quad v\_n^2 = \left(n - \frac{N}{2}\right) \Delta v^2, \quad m, n = 0, 1, \dots, N - 1, \dots$$

and

$$k\_p^1 = (p - \frac{N}{2})\Delta k^1, \quad k\_q^2 = (q - \frac{N}{2})\Delta k^2, \quad p, q = 0, 1, \dots, N - 1, \dots$$

where *N* is a power of 2, and the step sizes Δ*v*1, Δ*v*2, Δ*k*1, and Δ*k*<sup>2</sup> observe the Nyquist relations: Δ*v*1Δ*k*<sup>1</sup> = Δ*v*2Δ*k*<sup>2</sup> = 2*π*/*N*. Dempster & Hong (2002) shown that the numerical approximation now is given by

$$V\_T(k\_{p\prime}^1 k\_q^2) \approx \frac{e^{-\alpha\_1 k\_p^1 - \alpha\_2 k\_q^2}}{(2\pi)^2} \Psi\_T(k\_{p\prime}^1 k\_q^2)\sqrt{}$$

for *k*<sup>1</sup> *<sup>p</sup>* � *<sup>k</sup>*<sup>2</sup> *<sup>q</sup>*, where

12 Will-be-set-by-IN-TECH

**1**. The most direct extension of the Carr and Madan method to the multi assets model is to

*<sup>H</sup>*(*S*1*T*, *<sup>S</sup>*2*T*)=(*S*1*<sup>T</sup>* <sup>−</sup> *<sup>K</sup>*1)+(*S*2*<sup>T</sup>* <sup>−</sup> *<sup>K</sup>*2)+, where *K*<sup>1</sup> and *K*<sup>2</sup> are strike prices [see, e.g. Kwok *et al* (2010)]. Define *Xit* = log *Sit* and *ki* = log *Ki*, *i* = 1, 2, and let *fT*(*x*1, *x*2) be the joint density of *X*1*<sup>T</sup>* and *X*2*<sup>T</sup>* under the risk neutral probability. Then, the characteristic function of *X*1*<sup>T</sup>* and *X*2*<sup>T</sup>* is defined by the following

<sup>i</sup>(*u*<sup>1</sup> *<sup>x</sup>*1+*u*<sup>2</sup> *<sup>x</sup>*2) *fT*(*x*1, *x*2)*dx*1*dx*2.

*<sup>α</sup>*1*k*2+*α*2*k*2*VT*(*k*1, *k*2),

*p*+*v*2*k*<sup>2</sup> *q* )

*<sup>k</sup>*<sup>2</sup> )*fT*(*x*1, *x*2)*dx*1*dx*2.

.

*ψT*(*v*1, *v*2)*dv*1*dv*<sup>2</sup>

<sup>Δ</sup>*v*2, *<sup>m</sup>*, *<sup>n</sup>* <sup>=</sup> 0, 1, . . . , *<sup>N</sup>*−1,

<sup>Δ</sup>*k*2, *<sup>p</sup>*, *<sup>q</sup>* <sup>=</sup> 0, 1, . . . , *<sup>N</sup>*−1,

*<sup>p</sup>*, *<sup>k</sup>*<sup>2</sup> *q*), *<sup>n</sup>*)Δ*v*1Δ*v*2,

*<sup>k</sup>*<sup>1</sup> )(*ex*<sup>2</sup> <sup>−</sup> *<sup>e</sup>*

price the correlation option, whose payoff at the maturity *T* is defined by

 ∞ −∞  ∞ −∞ *e*

 ∞ *k*1

Following the Carr and Madam method we consider the modified call price:

*<sup>T</sup>*(*k*1, *k*<sup>2</sup> : *α*1, *α*2) = *e*

(*ex*<sup>1</sup> <sup>−</sup> *<sup>e</sup>*

for some parameters *α*<sup>1</sup> > 0 and *α*<sup>2</sup> > 0, which are chosen such that this modified price is square integrable for negative value of *k*<sup>1</sup> and *k*2. Then, by a direct integral the Fourier

*<sup>T</sup>*(*k*1, *k*<sup>2</sup> : *α*1, *α*2)

(*α*<sup>1</sup> + i*v*1)(*α*<sup>1</sup> + 1 + i*v*1)(*α*<sup>2</sup> + i*v*2)(*α*<sup>2</sup> + 1 + i*v*2)

Applying the Fourier inversion on *ψT*(*v*1, *v*2) and using the numerical approximation of

 ∞ −∞ *e* <sup>−</sup>i(*v*1*k*<sup>1</sup>

*N*−1 ∑ *n*=0 *e* <sup>−</sup>i(*v*<sup>1</sup> *mk*<sup>1</sup> *p*+*v*<sup>2</sup> *nk*2 *q* ) *ψT*(*v*<sup>1</sup> *<sup>m</sup>*, *<sup>v</sup>*<sup>2</sup>

*<sup>n</sup>* <sup>−</sup> *<sup>N</sup>* 2 

*<sup>q</sup>* <sup>−</sup> *<sup>N</sup>* 2 

where *N* is a power of 2, and the step sizes Δ*v*1, Δ*v*2, Δ*k*1, and Δ*k*<sup>2</sup> observe the Nyquist relations: Δ*v*1Δ*k*<sup>1</sup> = Δ*v*2Δ*k*<sup>2</sup> = 2*π*/*N*. Dempster & Hong (2002) shown that the numerical

> <sup>−</sup>*α*1*k*<sup>1</sup> *p*−*α*2*k*<sup>2</sup> *q* (2*π*)<sup>2</sup> <sup>Ψ</sup>*T*(*k*<sup>1</sup>

*v*<sup>1</sup> − i(*α*<sup>1</sup> + 1), *v*<sup>2</sup> − i(*α*<sup>2</sup> + 1)

**3.3 Fast Fourier transform methods - multi assets**

*ϕT*(*u*1, *u*2) =

*V*˜

−*rT* ∞ *k*2

*<sup>ψ</sup>T*(*v*1, *<sup>v</sup>*2) = <sup>F</sup>*V*˜

 ∞ −∞

*N*−1 ∑ *m*=0

for *p*, *q* = 0, 1, . . . , *N*−1. To apply the two-dimensional of the FFT algorithm we set

*<sup>n</sup>* =

*<sup>q</sup>* =

<sup>=</sup> *<sup>e</sup>*−*rT <sup>ϕ</sup><sup>T</sup>*

<sup>−</sup>*α*1*k*<sup>1</sup> *p*−*α*2*k*<sup>2</sup> *q*

(2*π*)<sup>2</sup>

(2*π*)<sup>2</sup>

Δ*v*1, *v*<sup>2</sup>

Δ*k*1, *k*<sup>2</sup>

*VT*(*k*<sup>1</sup> *<sup>p</sup>*, *<sup>k</sup>*<sup>2</sup> *<sup>q</sup>*) <sup>≈</sup> *<sup>e</sup>*

two-dimensional Fourier inversion integral we have

≈ *e* <sup>−</sup>*α*1*k*<sup>1</sup> *p*−*α*2*k*<sup>2</sup> *q*

*<sup>m</sup>* <sup>−</sup> *<sup>N</sup>* 2 

> *<sup>p</sup>* <sup>−</sup> *<sup>N</sup>* 2

two-dimensional Fourier transform:

transform of *V*˜(*k*1, *k*2) is given by

*VT*(*k*<sup>1</sup> *<sup>p</sup>*, *<sup>k</sup>*<sup>2</sup> *<sup>q</sup>*) = *<sup>e</sup>*

*v*1 *<sup>m</sup>* =

> *k*1 *<sup>p</sup>* =

approximation now is given by

and

Consider the price of the correlation option:

*VT*(*k*1, *k*2) = *e*

$$\Psi\_T(k\_{p\prime}^1 k\_q^2) = (-1)^{p+q} \sum\_{m=0}^{N-1} \sum\_{n=0}^{N-1} e^{-i\frac{2\pi}{N}(mp+nq)} \left( (-1)^{m+n} \psi\_T(v\_m^1, v\_n^2) \Delta v^1 \Delta v^2 \right)\_N$$

which can be computed via the FFT algorithm.

There are some other types of terminal payoff functions that admit analytic representation of the Fourier transform of the damped option price [see e.g. Eberlein *et al* (2010), Lee (2004), and references therein]. However, to derive the FFT pricing algorithm for the spread option:

$$H(\mathbf{S}\_{1T}, \mathbf{S}\_{2T}) = (\mathbf{S}\_{1T} - \mathbf{S}\_{2T} - \mathbf{K})^{+},\tag{13}$$

Dempster & Hong (2002) approximated the exercise region of *H*(*S*1*T*, *S*2*T*) by a combination of rectangular strips.

**2**. Hurd & Zhou (2010) proposed an alternative approach to pricing the European spread option (13) under the three-factor SV model and exponential Lévy models. Let *h*(*x*1, *x*2) be the terminal spread option payoff with the strike price *K* = 1, i.e.

$$h(\mathbf{x}\_1, \mathbf{x}\_2) = (e^{\mathbf{x}\_1} - e^{\mathbf{x}\_2} - 1)^+. \tag{14}$$

And denote Γ(*z*) the complex gamma function defined by

$$\Gamma(z) = \int\_0^\infty e^{-t} t^{z-1} dt, \quad \text{Re}\{z\} > 0.$$

Then, the method of Hurd and Zhou mainly relies on the following Fourier representation of the payoff function *h*(*x*1, *x*2):

$$h(\mathbf{x}\_1, \mathbf{x}\_2) = \frac{1}{(2\pi)^2} \int\_{\mathbf{i}\varepsilon\_2 - \infty}^{\mathbf{i}\varepsilon\_2 + \infty} \int\_{\mathbf{i}\varepsilon\_1 - \infty}^{\mathbf{i}\varepsilon\_1 + \infty} e^{\mathbf{i}(u\_1\mathbf{x}\_1 + u\_2\mathbf{x}\_2)} \hat{h}(u\_1, u\_2) du\_1 du\_2. \tag{15}$$

where *�*<sup>1</sup> and *�*<sup>2</sup> are any two given real numbers with *�*<sup>2</sup> > 0 and *�*<sup>1</sup> + *�*<sup>2</sup> < −1, and

$$
\hat{h}(\boldsymbol{\mu}\_1, \boldsymbol{\mu}\_2) = \frac{\Gamma\left(\mathrm{i}(\boldsymbol{\mu}\_1 + \boldsymbol{\mu}\_2) - 1\right)\Gamma(-\mathrm{i}\boldsymbol{\mu}\_2)}{\Gamma(\mathrm{i}\boldsymbol{\mu}\_1 + 1)}.
$$

The proof of this formula can be found in the appendix of Hurd & Zhou (2010). Let *X*1*<sup>t</sup>* = log *S*1*<sup>t</sup>* and *X*2*<sup>t</sup>* = log *S*2*<sup>t</sup>* with *X*<sup>10</sup> = *x*<sup>1</sup> and *X*<sup>20</sup> = *x*2, and let *ϕT*(*u*1, *u*2) be the characteristic function of *X*1*<sup>T</sup>* − *x*<sup>1</sup> and *X*2*<sup>T</sup>* − *x*2. Then, we have

$$\mathbb{E}\left[e^{\mathbf{i}\left(\mu\_1 X\_{17} + \mu\_2 X\_{27}\right)}\right] = e^{\mathbf{i}\left(\mu\_1 \mathbf{x}\_1 + \mu\_2 \mathbf{x}\_2\right)} \mathbf{q}\_T(\mu\_1, \mu\_2).$$

Now, using the formula (16), the price of the spread option (14) is expressed as an explicit two-dimensional Fourier inversion transform:

$$\begin{split} V\_{T}(\mathbf{x}\_{1},\mathbf{x}\_{2}) &= e^{-rT} \mathbb{E} \left[ \left( e^{\mathbf{x}\_{1\mathrm{T}}} - e^{\mathbf{x}\_{2\mathrm{T}}} - 1 \right)^{+} \right] \\ &= e^{-rT} \mathbb{E} \left[ \frac{1}{(2\pi)^{2}} \int\_{i\mathbf{i}\varepsilon\_{2} - \infty}^{i\varepsilon\_{2} + \infty} \int\_{i\mathbf{i}\varepsilon\_{1} - \infty}^{i\varepsilon\_{1} + \infty} e^{i(\mu\_{1}\mathbf{X}\_{1\mathrm{T}} + \mu\_{2}\mathbf{X}\_{2\mathrm{T}})} \hat{h}(\mu\_{1}, \mu\_{2}) du\_{1} du\_{2} \right] \\ &= \frac{e^{-rT}}{(2\pi)^{2}} \int\_{i\mathbf{i}\varepsilon\_{2} - \infty}^{i\varepsilon\_{2} + \infty} \int\_{i\mathbf{i}\varepsilon\_{1} - \infty}^{i\varepsilon\_{1} + \infty} \mathbb{E} \left[ e^{i(\mu\_{1}\mathbf{X}\_{1\mathrm{T}} + \mu\_{2}\mathbf{X}\_{2\mathrm{T}})} \right] \hat{h}(\mu\_{1}, \mu\_{2}) du\_{1} du\_{2} \\ &= \frac{e^{-rT}}{(2\pi)^{2}} \int\_{i\mathbf{i}\varepsilon\_{2} - \infty}^{i\varepsilon\_{2} + \infty} \int\_{i\mathbf{i}\varepsilon\_{1} - \infty}^{i\varepsilon\_{1} + \infty} e^{i(\mu\_{1}\mathbf{x}\_{1} + \mu\_{2}\mathbf{x}\_{2})} \, q\_{T}(\mu\_{1}, \mu\_{2}) \hat{h}(\mu\_{1}, \mu\_{2}) du\_{1} du\_{2} . \end{split}$$

where Φ*n*(0, *b*) and Ψ*n*(0, *b*) are two integrals given by:

*<sup>e</sup><sup>x</sup>* cos �

for any [*c*, *d*] ⊂ [*a*, *b*], which are analytically given by

+ *nπ b* − *a*

⎧ ⎨ ⎩

integration [*a*, *b*] in that paper, which is given by

[*a*, *<sup>b</sup>*] = �

� sin � *nπ d* − *a b* − *a*

1 + � *<sup>n</sup><sup>π</sup> b*−*a* �2 � cos � *nπ d* − *a b* − *a*

> sin � *nπ d* − *a b* − *a*

*c*<sup>1</sup> − *δ* �

numerical comparisons for these different methods is also given in their paper.

**4. Fourier transform methods for pricing path dependent options**

monitored dates may be many times more than the exercise dates. Denote

where *K* is the strike price, *S* is the spot price of the underlying asset. Let

*tmL* : *m* = 1, . . . , *M*−1

**4.1 The CONV method for pricing Bermudan barrier options**

price Bermudan options under exponential Lévy models.

� **T** = �

**T***<sup>e</sup>* = �

*nπ x* − *a b* − *a* �

*dx* and <sup>Ψ</sup>*n*(*c*, *<sup>d</sup>*) = � *<sup>d</sup>*

� *e*

*nπ c* − *a b* − *a*

(*d* − *c*) *n* = 0

� *e <sup>d</sup>* <sup>−</sup> *<sup>n</sup><sup>π</sup> b* − *a*

Fourier Transform Methods for Option Pricing 279

<sup>−</sup> sin �

Fang & Oosterlee (2008) also showed that, in most cases, the convergence rate of the COS formula (20) is exponential and the computational complexity is linear. They also discussed the truncation range for COS method, and gave a general formula to determine the interval of

*<sup>c</sup>*<sup>2</sup> <sup>+</sup> <sup>√</sup>*c*4, *<sup>c</sup>*<sup>1</sup> <sup>+</sup> *<sup>δ</sup>*

where *c*1, *c*2, and *c*<sup>4</sup> are the first, second, and fourth cumulates of *Xt* = log(*St*/*K*). Also, the constant *δ* depends on the tolerance level in the approximation (17), and usually we choose *δ* = 10. Meanwhile, Ding & U (2011) respectively applied the Fourier-sine and Fourier expansions to substitute the Fourier-cosine expansion in (18). A comprehensive analysis with

Pricing Bermudan or barrier options is much harder than pricing European options. Because these options are depended on paths of the price process for the underlying assets. Recently, some new numerical integration methods based on Fourier transforms are proposed. Lord *et al* (2008) proposed an efficient and accurate FFT-based method, called the CONV method, to

In the following we apply the CONV method to price a Bermudan barrier option in which the

*<sup>G</sup>*(*S*) = � (*<sup>S</sup>* <sup>−</sup> *<sup>K</sup>*)+, for call option,

(*<sup>K</sup>* <sup>−</sup> *<sup>S</sup>*)+, for put option,

*tmL*<sup>+</sup>*<sup>l</sup>* : *m* = 0, 1, . . . , *M*, *l* = 0, 1, . . . , *L*−1

�

*<sup>d</sup>* <sup>−</sup> cos �

�

*<sup>c</sup>*<sup>2</sup> <sup>+</sup> <sup>√</sup>*c*<sup>4</sup>

�

� ,

� <sup>⊂</sup> **<sup>T</sup>**, (22)

sin � *nπ c* − *a b* − *a*

*c*

*nπ c* − *a b* − *a*

�� *b* − *a nπ*

cos � *nπ x* − *a b* − *a*

> � *e c*

� *e c* � ,

*n* � 0,

.

, (21)

� *dx*,

<sup>Φ</sup>*n*(*c*, *<sup>d</sup>*) = � *<sup>d</sup>*

*c*

<sup>Φ</sup>*n*(*c*, *<sup>d</sup>*) = <sup>1</sup>

Ψ*n*(*c*, *d*) =

This two-dimensional Fourier inversion integral can be numerically computed by the usual FFT calculations. Since this approach does not require an approximation of the exercise region, it is considered to be more computationally efficient.

#### **3.4 Fourier expansion methods**

Fang & Oosterlee (2008) proposed a new numerical method for pricing European options, which is called the COS method. This method is based on the Fourier transform and the Fourier-cosine expansion, and is shown to have the exponential convergence rate and linear computational complexity. And now, this novel method has been considered by many authors to price various options [see e.g. Ding *et al* (2011b),Fang & Oosterlee (2009), and the references therein].

We recall that the price *St* of underlying asset follows an exponential Lévy model. Let *Xt* = log(*St*/*K*), where *K* is the strike price. Consider the price of European call in the form:

$$V\_0(K) = e^{-rT} \mathbb{E}\left[ (\mathbb{S}\_T - K)^+ \right] = e^{-rT} K \int\_0^\infty (e^\chi - 1) f\_T(\mathbf{x}) d\mathbf{x} \tag{16}$$

where *fT*(*x*) is the probability density of *XT* under the risk-neutral probability. Let *ϕT*(*u*) be the characteristic function of *fT*(*x*). The main ideas of the COS method are to choose two numbers *a* and *b* such that the truncated integral approximates the infinite integral very well, i.e.,

$$\tilde{\varphi}\_{\mathsf{T}}(\mathsf{u}) = \int\_{a}^{b} \mathsf{e}^{\mathsf{i}\mathsf{u}\mathsf{x}} f\_{\mathsf{T}}(\mathsf{x}) d\mathsf{x} \approx \int\_{-\infty}^{\infty} \mathsf{e}^{\mathsf{i}\mathsf{u}\mathsf{x}} f\_{\mathsf{T}}(\mathsf{x}) d\mathsf{x} = \varphi\_{\mathsf{T}}(\mathsf{u}),\tag{17}$$

and to consider the Fourier-cosine expansion of *fT*(*x*) in [*a*, *b*]:

$$f\_T(\mathbf{x}) = \frac{1}{2}A\_0 + \sum\_{n=1}^{\infty} A\_n \cos \left( n\pi \frac{\mathbf{x} - a}{b - a} \right), \quad \mathbf{x} \in [a, b], \tag{18}$$

where

$$A\_n = \frac{2}{b-a} \int\_a^b \mathfrak{p}\_T(\mathbf{x}) \cos \left( n\pi \frac{\mathbf{x} - a}{b-a} \right) d\mathbf{x}, \quad n = 0, 1, 2, \dots$$

From (17) we have

$$A\_{\rm H} = \frac{2}{b-a} \text{Re} \{ \vec{\phi}\_T \left( \frac{n\pi}{b-a} \right) \exp \left( -\mathbf{i} \frac{na\pi}{b-a} \right) \}, \quad n = 0, 1, 2, \dots, 4$$

Then, we get an approximation of *fT*(*x*) in (18) by

$$f\_T(\mathbf{x}) \approx \frac{1}{2}\tilde{A}\_0 + \sum\_{n=1}^{N-1} \tilde{A}\_n \cos\left(n\pi \frac{\mathbf{x} - a}{b - a}\right), \quad \mathbf{x} \in [a, b]. \tag{19}$$

where

$$\tilde{A}\_n = \frac{2}{b-a} \text{Re} \left\{ \varphi\_T \left( \frac{n\pi}{b-a} \right) \exp \left( -\mathbf{i} \frac{na\pi}{b-a} \right) \right\} \approx A\_{\text{II}} \quad n = 0, 1, \ldots, N - 1$$

Now, substituting (19) into (16), we obtain an approximation of the option price:

$$V\_0(K) \approx K e^{-rT} \left\{ \frac{1}{2} \tilde{A}\_0 \left( \Phi\_0(0, b) - \Psi\_0(0, b) \right) + \sum\_{n=1}^{N-1} \tilde{A}\_n \left( \Phi\_n(0, b) - \Psi\_n(0, b) \right) \right\},\tag{20}$$

14 Will-be-set-by-IN-TECH

This two-dimensional Fourier inversion integral can be numerically computed by the usual FFT calculations. Since this approach does not require an approximation of the exercise region,

Fang & Oosterlee (2008) proposed a new numerical method for pricing European options, which is called the COS method. This method is based on the Fourier transform and the Fourier-cosine expansion, and is shown to have the exponential convergence rate and linear computational complexity. And now, this novel method has been considered by many authors to price various options [see e.g. Ding *et al* (2011b),Fang & Oosterlee (2009), and the references

We recall that the price *St* of underlying asset follows an exponential Lévy model. Let *Xt* = log(*St*/*K*), where *K* is the strike price. Consider the price of European call in the form:

> = *e* <sup>−</sup>*rTK*

> > ∞ −∞ *e*

where *fT*(*x*) is the probability density of *XT* under the risk-neutral probability. Let *ϕT*(*u*) be the characteristic function of *fT*(*x*). The main ideas of the COS method are to choose two numbers *a* and *b* such that the truncated integral approximates the infinite integral very well,

 ∞ 0

 + *N*−1 ∑ *n*=1 *A*˜*n* 

d*x*, *n* = 0, 1, 2, . . . .

, *n* = 0, 1, 2, . . . .

≈ *An*, *n* = 0, 1, . . . , *N*−1

Φ*n*(0, *b*) − Ψ*n*(0, *b*)

(*e<sup>x</sup>* <sup>−</sup> <sup>1</sup>)*fT*(*x*)*dx*, (16)

<sup>i</sup>*ux fT*(*x*)*dx* = *ϕT*(*u*), (17)

, *x* ∈ [*a*, *b*], (18)

, *x* ∈ [*a*, *b*], (19)

, (20)

(*ST* <sup>−</sup> *<sup>K</sup>*)<sup>+</sup>

<sup>i</sup>*ux fT*(*x*)*dx* <sup>≈</sup>

∞ ∑ *n*=1

*p*˜*T*(*x*) cos

*N*−1 ∑ *n*=1

 exp − i *naπ b* − *a*

Now, substituting (19) into (16), we obtain an approximation of the option price:

Φ0(0, *b*) − Ψ0(0, *b*)

*An* cos *nπ x* − *a b* − *a*

> *nπ x* − *a b* − *a*

*A*˜*<sup>n</sup>* cos *nπ x* − *a b* − *a*

 exp − i *naπ b* − *a*

it is considered to be more computationally efficient.

*V*0(*K*) = *e*

*ϕ*˜*T*(*u*) =

*fT*(*x*) = <sup>1</sup>

*b* − *a*

*An* <sup>=</sup> <sup>2</sup>

*b* − *a* Re *ϕ*˜*T nπ b* − *a*

Then, we get an approximation of *fT*(*x*) in (18) by

*fT*(*x*) <sup>≈</sup> <sup>1</sup>

2 *A*˜ 0 

*An* <sup>=</sup> <sup>2</sup>

*<sup>A</sup>*˜*<sup>n</sup>* <sup>=</sup> <sup>2</sup>

*b* − *a* Re *ϕT nπ b* − *a*

*<sup>V</sup>*0(*K*) <sup>≈</sup> *Ke*−*rT*<sup>1</sup>

<sup>−</sup>*rT***E** 

> *b a e*

and to consider the Fourier-cosine expansion of *fT*(*x*) in [*a*, *b*]:

<sup>2</sup> *<sup>A</sup>*<sup>0</sup> <sup>+</sup>

 *b a*

<sup>2</sup> *<sup>A</sup>*˜ <sup>0</sup> <sup>+</sup>

**3.4 Fourier expansion methods**

therein].

i.e.,

where

where

From (17) we have

where Φ*n*(0, *b*) and Ψ*n*(0, *b*) are two integrals given by:

$$\Phi\_n(c,d) = \int\_{\mathcal{c}}^d e^{\chi} \cos \left(n\pi \frac{\chi-a}{b-a}\right) d\chi \quad \text{and} \quad \Psi\_n(c,d) = \int\_{\mathcal{c}}^d \cos \left(n\pi \frac{\chi-a}{b-a}\right) d\chi$$

for any [*c*, *d*] ⊂ [*a*, *b*], which are analytically given by

$$\Phi\_{\rm ll}(c,d) = \frac{1}{1 + \left(\frac{n\pi}{b-a}\right)^2} \left[ \cos\left(n\pi\frac{d-a}{b-a}\right)e^d - \cos\left(n\pi\frac{c-a}{b-a}\right)e^c \right]$$

$$+ \frac{n\pi}{b-a} \sin\left(n\pi\frac{d-a}{b-a}\right)e^d - \frac{n\pi}{b-a} \sin\left(n\pi\frac{c-a}{b-a}\right)e^c \Big|\_{\prime}$$

$$\Psi\_{\rm n}(c,d) = \begin{cases} \left[ \sin\left(n\pi\frac{d-a}{b-a}\right) - \sin\left(n\pi\frac{c-a}{b-a}\right) \right] \frac{b-a}{n\pi} & n \neq 0, \\\ (d-c) & n=0 \end{cases}$$

Fang & Oosterlee (2008) also showed that, in most cases, the convergence rate of the COS formula (20) is exponential and the computational complexity is linear. They also discussed the truncation range for COS method, and gave a general formula to determine the interval of integration [*a*, *b*] in that paper, which is given by

$$\mathbf{a}[a,b] = \left[c\_1 - \delta\sqrt{c\_2 + \sqrt{c\_4}}, c\_1 + \delta\sqrt{c\_2 + \sqrt{c\_4}}\right],\tag{21}$$

where *c*1, *c*2, and *c*<sup>4</sup> are the first, second, and fourth cumulates of *Xt* = log(*St*/*K*). Also, the constant *δ* depends on the tolerance level in the approximation (17), and usually we choose *δ* = 10. Meanwhile, Ding & U (2011) respectively applied the Fourier-sine and Fourier expansions to substitute the Fourier-cosine expansion in (18). A comprehensive analysis with numerical comparisons for these different methods is also given in their paper.

#### **4. Fourier transform methods for pricing path dependent options**

#### **4.1 The CONV method for pricing Bermudan barrier options**

Pricing Bermudan or barrier options is much harder than pricing European options. Because these options are depended on paths of the price process for the underlying assets. Recently, some new numerical integration methods based on Fourier transforms are proposed. Lord *et al* (2008) proposed an efficient and accurate FFT-based method, called the CONV method, to price Bermudan options under exponential Lévy models.

In the following we apply the CONV method to price a Bermudan barrier option in which the monitored dates may be many times more than the exercise dates. Denote

$$G(\mathcal{S}) = \begin{cases} (\mathcal{S} - K)^+ \text{, for call option,} \\ (K - \mathcal{S})^+ \text{, for put option,} \end{cases}$$

where *K* is the strike price, *S* is the spot price of the underlying asset. Let

$$\begin{cases} \mathbb{T}\_{\cdot} = \{t\_{mL+l} : m = 0, 1, \dots, M, \ l = 0, 1, \dots, L-1\}, \\\mathbb{T}\_{\ell} = \{t\_{mL} : m = 1, \dots, M-1\} \subset \mathbb{T}\_{\prime} \end{cases} \tag{22}$$

where *<sup>v</sup>*(*x*, *tk*) = *<sup>V</sup>*(*S*0*ex*, *tk*) for any *tk* <sup>∈</sup> **<sup>T</sup>**, and *<sup>h</sup>* <sup>=</sup> log(*H*/*S*0).

*f*(· | *x*) possesses the property:

*c*(*x*, *tk*) = *e*

(Convolution), the integral becomes

*e r*Δ*t* F *e*

a simple calculation, we have

*e*

Denote

in (26):

*e r*Δ*t* F *e*

*<sup>α</sup>xc*(*x*, *tk*) = *e*

function *f*(−*x*), i.e.

Thus we have

infinite integrals *c*(*x*, *tk*) in (24) it becomes to

−*r*Δ*t* ∞ −∞

*e*

Since each Lévy process is stationary and has independent increments, the condition density

Fourier Transform Methods for Option Pricing 281

*f*(*y*|*x*) = *f*(*y*−*x*), *x*, *y* ∈ **R**,

where *f*(*y*) is the density of *Xt*<sup>1</sup> under the initial condition *Xt*<sup>0</sup> = *x*. Applying this property to

for any *tk* ∈ **T**. Then, this integral can be rewritten as a convolution of *v*(*x*+*z*, *tk*<sup>+</sup>1) and the

where *α* > 0 is chosen to improve the integrability. Now, the main idea of CONV method is that, taking the Fourier transform in the both sides of (25) and applying the Property P3

Denote *ϕ*(*z*) is the characteristic function of the density *f*(*x*) on the complex plan **C**. Then, by

−(*u*−i*α*)

−(*u*−i*α*)

*<sup>α</sup><sup>x</sup> <sup>f</sup>*(−*x*)

<sup>−</sup>*r*Δ*<sup>t</sup> <sup>f</sup>*(−*x*) <sup>∗</sup> *<sup>v</sup>*(*x*, *tk*<sup>+</sup>1).

 ∗ *e*

> ∗ *e*

 · F *e*

> <sup>∞</sup> −∞ *e*

 *∂ϕT*(*u*) *∂u*

Δ*x*, *j* = 0, 1, . . . , *N*−1,

*N*)Δ*u*, *yj* = *xj*, *j* = 0, 1, . . . , *N*−1,

 *u*=0 2 ,

 (*u*) · F *e*

−*r*Δ*t* ∞ −∞

*<sup>α</sup>xv*(*x*, *tk*<sup>+</sup>1)

*<sup>α</sup>xv*(*x*, *tk*<sup>+</sup>1)

*<sup>α</sup>xv*(*x*, *tk*<sup>+</sup>1)

*<sup>α</sup>xv*(*x*, *tk*<sup>+</sup>1)

(*u*)

 (*u*).

<sup>i</sup>*uy*+*αyv*(*y*, *tk*<sup>+</sup>1)*dydu*. (26)

 (*u*).

*v*(*x*+*z*, *tk*<sup>+</sup>1)*f*(*z*)*dz*

, (25)

*v*(*y*, *tk*<sup>+</sup>1)*f*(*y*−*x*)*dy* = *e*

*c*(*x*, *tk*) = *e*

−*r*Δ*t e <sup>α</sup><sup>x</sup> <sup>f</sup>*(−*x*)

(*u*) = <sup>F</sup>*<sup>e</sup>*

<sup>=</sup> <sup>F</sup> *e <sup>α</sup><sup>x</sup> <sup>f</sup>*(−*x*)

*<sup>α</sup>xc*(*x*, *tk*) = *e*

*<sup>α</sup>xc*(*x*, *tk*)

*<sup>α</sup>xc*(*x*, *tk*)

Thus, taking the inverse Fourier transform we obtain

<sup>−</sup>*r*Δ*<sup>t</sup>* <sup>1</sup> 2*π*

*κ* = *δ*

*xj* =

2

*uj* = (*<sup>j</sup>* <sup>−</sup> <sup>1</sup>

*<sup>j</sup>* <sup>−</sup> <sup>1</sup> 2 *N*

 (*u*) = *ϕ* 

 ∞ −∞ *e* <sup>−</sup>i*uxϕ* 

−*∂*2*ϕT*(*u*) *∂u*<sup>2</sup>

 *u*=0 +

where *ϕT*(*u*) is the characteristic function of *XT*, and *δ* is a proportionality constant. According to the suggestion from Lord *et al* (2008), we can take *δ* = 20 for the GBM model and *δ* = 40 for other exponential Lévy models. Let *N* be a power of 2. We consider the grid points on *x*-axes:

where Δ*x* = *κ*/*N*. Furthermore, we also consider the grid points for the numerical integrals

be the set of pre-specified monitored dates and the set of pre-specified exercise dates, respectively, before maturity *T*, where

$$0 = t\_0 < t\_1 < \dots < t\_{ML} = T \quad \text{with} \quad \Delta t = t\_k - t\_{k-1} = T/(ML).$$

Consider a discrete American barrier option, which is monitored at every *tk* ∈ **T** and can be exercised at each *tk* ∈ **T***e*, namely Bermudan barrier option, whose payoff is given by

$$\mathcal{G}(\mathcal{S}\_{l\_k})\mathbf{1}\_{\{\mathcal{S}\_{l\_k} < H\}} + \mathcal{R}\_0 \mathbf{1}\_{\{\mathcal{S}\_{l\_k} \ge H\}}\prime \quad t\_k \in \mathbb{T}.$$

Here *Stk* is the price of the underlying asset at time *tk* ∈ **T**, *H* > *K* is the constant barrier and *R*<sup>0</sup> is the contractual rebate. That is, this Bermudan barrier option is an up-and-out barrier option that cease to exist if the asset price *Stk* hits the barrier level *H* at one time *tk* ∈ **T**, and it can also be exercised at any time *tk* ∈ **T***<sup>e</sup>* before maturity *T*.

Denote *V*(*S*, *tk*) the value of this Bermudan barrier option at time *tk* and the spot price *Stk* = *S*. With help of the risk-neutral valuation formula, this price process can be computed recursively by the following backward induction:

$$\begin{cases} V(\mathcal{S}, t\_{ML}) = G(\mathcal{S}) \mathbf{1}\_{\{\mathcal{S} < H\}} + \mathcal{R}\_0 \mathbf{1}\_{\{\mathcal{S} \ge H\}}, \\ \mathcal{C}(\mathcal{S}, t\_k) \\ V(\mathcal{S}, t\_k) = \mathcal{C}(\mathcal{S}, t\_k) \mathbf{1}\_{\{\mathcal{S} < H\}} + e^{-r(T - t\_k)} \mathcal{R}\_0 \mathbf{1}\_{\{\mathcal{S} \ge H\}}, \quad t\_k \in \mathbb{T} \backslash \mathsf{T}\_{\mathcal{C}} \\ V(\mathcal{S}, t\_k) = \max\left\{ G(\mathcal{S}), \mathcal{C}(\mathcal{S}, t\_k) \right\} \mathbf{1}\_{\{\mathcal{S} < H\}} + e^{-r(T - t\_k)} \mathcal{R}\_0 \mathbf{1}\_{\{\mathcal{S} \ge H\}}, \quad t\_k \in \mathsf{T}\_{\mathcal{C}} \end{cases} \tag{23}$$

in specialty, the initial price is given

$$V(\mathcal{S}, t\_0) = \mathcal{C}(\mathcal{S}, t\_0) = \mathbb{E}\left[e^{-r\Delta t} V(\mathcal{S}, t\_1) \mid \mathcal{S}\_{t\_0} = \mathcal{S}\right].$$

Here *r* > 0 is the interest rate, and **E**[·|·] is the conditional expectation under the risk-neutral probability.

Assume that the price process of the underlying asset is given by

$$S\_t = S\_0 e^{X\_t}, \quad t \ge 0,$$

where *Xt* is a Lévy process and *S*<sup>0</sup> is the initial price. Let *f*(· | *x*) be the condition density of *Xtk*<sup>+</sup><sup>1</sup> given *Xtk* = *x* for *tk* ∈ **T**. Set

$$g(\boldsymbol{x}) = \begin{cases} (\mathbb{S}\_0 \boldsymbol{\varepsilon}^{\boldsymbol{\chi}} - \boldsymbol{K})^+, & \text{for a call option,} \\ (\boldsymbol{K} - \mathbb{S}\_0 \boldsymbol{\varepsilon}^{\boldsymbol{\chi}})^+, & \text{for a put option.} \end{cases}$$

Then the backward induction (23) can be rewritten by

$$\begin{cases} v(\mathbf{x}, t\_{ML}) = g(\mathbf{x}) \mathbf{1}\_{\{\mathbf{x} < h\}} + R\_0 \mathbf{1}\_{\{\mathbf{x} \ge h\}}, \\ c(\mathbf{x}, t\_k) \\ \begin{cases} v(\mathbf{x}, t\_k) \\ v(\mathbf{x}, t\_k) \end{cases} = e(\mathbf{x}, t\_k) \mathbf{1}\_{\{\mathbf{x} < h\}} + R\_0 \mathbf{1}\_{\{\mathbf{x} \ge h\}}, \quad t\_k \in \mathbb{T} \backslash \mathbb{T}\_{\mathcal{E}\_\ell} \\ v(\mathbf{x}, t\_k) \quad = \max\left\{ g(\mathbf{x}), c(\mathbf{x}, t\_k) \right\} \mathbf{1}\_{\{\mathbf{x} < h\}} + e^{-r(T - t\_k)} R\_0 \mathbf{1}\_{\{\mathbf{x} \ge h\}}, \quad t\_k \in \mathbb{T}\_{\mathcal{E}\_\ell} \end{cases} \tag{24}$$

and

$$v(\mathbf{x}, t\_0) = c(\mathbf{x}, t\_0) = e^{-r\Delta t} \int\_{-\infty}^{\infty} v(y, t\_1) f(y \mid \mathbf{x}) \, \mathbf{d}y \, \mathbf{x}$$

where *<sup>v</sup>*(*x*, *tk*) = *<sup>V</sup>*(*S*0*ex*, *tk*) for any *tk* <sup>∈</sup> **<sup>T</sup>**, and *<sup>h</sup>* <sup>=</sup> log(*H*/*S*0).

Since each Lévy process is stationary and has independent increments, the condition density *f*(· | *x*) possesses the property:

$$f(y|\mathbf{x}) = f(y-\mathbf{x}), \quad \mathbf{x}, y \in \mathbb{R}\_{\succ}$$

where *f*(*y*) is the density of *Xt*<sup>1</sup> under the initial condition *Xt*<sup>0</sup> = *x*. Applying this property to infinite integrals *c*(*x*, *tk*) in (24) it becomes to

$$\mathcal{L}(\mathbf{x}, t\_k) = e^{-r\Delta t} \int\_{-\infty}^{\infty} v(y, t\_{k+1}) f(y-\mathbf{x}) dy = e^{-r\Delta t} \int\_{-\infty}^{\infty} v(\mathbf{x} + \mathbf{z}, t\_{k+1}) f(\mathbf{z}) d\mathbf{z}$$

for any *tk* ∈ **T**. Then, this integral can be rewritten as a convolution of *v*(*x*+*z*, *tk*<sup>+</sup>1) and the function *f*(−*x*), i.e.

$$
\sigma(\mathfrak{x}, t\_k) = e^{-r\Delta t} f(-\mathfrak{x}) \* \upsilon(\mathfrak{x}, t\_{k+1}) .
$$

Thus we have

16 Will-be-set-by-IN-TECH

be the set of pre-specified monitored dates and the set of pre-specified exercise dates,

<sup>0</sup> = *<sup>t</sup>*<sup>0</sup> < *<sup>t</sup>*<sup>1</sup> < ··· < *tML* = *<sup>T</sup>* with <sup>Δ</sup>*<sup>t</sup>* = *tk* − *tk*−<sup>1</sup> = *<sup>T</sup>*/(*ML*). Consider a discrete American barrier option, which is monitored at every *tk* ∈ **T** and can be

Here *Stk* is the price of the underlying asset at time *tk* ∈ **T**, *H* > *K* is the constant barrier and *R*<sup>0</sup> is the contractual rebate. That is, this Bermudan barrier option is an up-and-out barrier option that cease to exist if the asset price *Stk* hits the barrier level *H* at one time *tk* ∈ **T**, and it

Denote *V*(*S*, *tk*) the value of this Bermudan barrier option at time *tk* and the spot price *Stk* = *S*. With help of the risk-neutral valuation formula, this price process can be computed

*V*(*Stk*<sup>+</sup><sup>1</sup> , *tk*<sup>+</sup>1) | *Stk* = *S*

*<sup>V</sup>*(*S*, *tk*) = *<sup>C</sup>*(*S*, *tk*)1{*S*<*H*} <sup>+</sup> *<sup>e</sup>*−*r*(*T*−*tk* )*R*01{*S*≥*H*}, *tk* <sup>∈</sup> **<sup>T</sup>** \ **<sup>T</sup>***e*,

�

� *e* −*r*Δ*t*

Here *r* > 0 is the interest rate, and **E**[·|·] is the conditional expectation under the risk-neutral

*St* <sup>=</sup> *<sup>S</sup>*0*eXt* , *<sup>t</sup>* <sup>≥</sup> 0,

where *Xt* is a Lévy process and *S*<sup>0</sup> is the initial price. Let *f*(· | *x*) be the condition density of

� (*S*0*e<sup>x</sup>* <sup>−</sup> *<sup>K</sup>*)+, for a call option, (*<sup>K</sup>* <sup>−</sup> *<sup>S</sup>*0*ex*)+, for a put option.

<sup>−</sup><sup>∞</sup> *<sup>v</sup>*(*y*, *tk*<sup>+</sup>1)*f*(*<sup>y</sup>* <sup>|</sup> *<sup>x</sup>*)d*y*, *tk* <sup>∈</sup> **<sup>T</sup>**,

�

−*r*Δ*t* � ∞ −∞

*G*(*S*), *C*(*S*, *tk*)

*V*(*S*, *t*0) = *C*(*S*, *t*0) = **E**

Assume that the price process of the underlying asset is given by

*g*(*x*) =

*<sup>v</sup>*(*x*, *tML*) = *<sup>g</sup>*(*x*)1{*x*<*h*} <sup>+</sup> *<sup>R</sup>*01{*x*≥*h*},

*<sup>v</sup>*(*x*, *tk*) = *<sup>c</sup>*(*x*, *tk*)1{*x*<*h*} <sup>+</sup> *<sup>R</sup>*01{*x*≥*h*}, *tk* <sup>∈</sup> **<sup>T</sup>** \ **<sup>T</sup>***e*,

*g*(*x*), *c*(*x*, *tk*)

*v*(*x*, *t*0) = *c*(*x*, *t*0) = *e*

Then the backward induction (23) can be rewritten by

*<sup>c</sup>*(*x*, *tk*) = *<sup>e</sup>*−*r*Δ*<sup>t</sup>* � <sup>∞</sup>

*<sup>v</sup>*(*x*, *tk*) = max �

*<sup>k</sup>* <sup>≥</sup>*H*}, *tk* <sup>∈</sup> **<sup>T</sup>**.

�

, *tk* ∈ **T**,

*V*(*S*, *t*1) | *St*<sup>0</sup> = *S*

<sup>1</sup>{*S*<*H*} <sup>+</sup> *<sup>e</sup>*−*r*(*T*−*tk* )*R*01{*S*≥*H*}, *tk* <sup>∈</sup> **<sup>T</sup>***e*,

<sup>1</sup>{*x*<*h*} <sup>+</sup> *<sup>e</sup>*−*r*(*T*−*tk* )*R*01{*x*≥*h*}, *tk* <sup>∈</sup> **<sup>T</sup>***e*,

*v*(*y*, *t*1)*f*(*y* | *x*)d*y*,

� . (23)

(24)

exercised at each *tk* ∈ **T***e*, namely Bermudan barrier option, whose payoff is given by

*<sup>k</sup>* <sup>&</sup>lt;*H*} <sup>+</sup> *<sup>R</sup>*01{*St*

*<sup>G</sup>*(*Stk* )1{*St*

can also be exercised at any time *tk* ∈ **T***<sup>e</sup>* before maturity *T*.

*<sup>V</sup>*(*S*, *tML*) = *<sup>G</sup>*(*S*)1{*S*<*H*} <sup>+</sup> *<sup>R</sup>*01{*S*≥*H*},

recursively by the following backward induction:

� *e*−*r*Δ*<sup>t</sup>*

*C*(*S*, *tk*) = **E**

in specialty, the initial price is given

*Xtk*<sup>+</sup><sup>1</sup> given *Xtk* = *x* for *tk* ∈ **T**. Set

⎧ ⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎩

and

*<sup>V</sup>*(*S*, *tk*) = max �

⎧ ⎪⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎪⎩

probability.

respectively, before maturity *T*, where

$$e^{a\mathbf{x}}c(\mathbf{x},t\_k) = e^{-r\Delta t} \left[ e^{a\mathbf{x}}f(-\mathbf{x}) \right] \* \left[ e^{a\mathbf{x}}v(\mathbf{x},t\_{k+1}) \right]\_{\prime} \tag{25}$$

where *α* > 0 is chosen to improve the integrability. Now, the main idea of CONV method is that, taking the Fourier transform in the both sides of (25) and applying the Property P3 (Convolution), the integral becomes

$$\begin{split} \left[e^{r\Delta t}\mathcal{F}\left[e^{\alpha x}c(\mathbf{x},t\_{k})\right](\boldsymbol{u}) &= \mathcal{F}\left\{\left[e^{\alpha x}f(-\mathbf{x})\right]\*\left[e^{\alpha x}v(\mathbf{x},t\_{k+1})\right]\right\}(\boldsymbol{u})\\ &= \mathcal{F}\left[e^{\alpha x}f(-\mathbf{x})\right](\boldsymbol{u})\cdot\mathcal{F}\left[e^{\alpha x}v(\mathbf{x},t\_{k+1})\right](\boldsymbol{u}). \end{split}$$

Denote *ϕ*(*z*) is the characteristic function of the density *f*(*x*) on the complex plan **C**. Then, by a simple calculation, we have

$$e^{r\Delta t} \mathcal{F}\left[e^{a\mathbf{x}}c(\mathbf{x},t\_k)\right](\mathsf{u}) = \varphi\left(-(\mathsf{u}-\mathsf{i}\mathsf{u})\right)\cdot \mathcal{F}\left[e^{a\mathbf{x}}v(\mathbf{x},t\_{k+1})\right](\mathsf{u}).$$

Thus, taking the inverse Fourier transform we obtain

$$e^{\mathrm{ax}}c(\mathbf{x},t\_k) = e^{-r\Delta t} \frac{1}{2\pi} \int\_{-\infty}^{\infty} e^{-\mathrm{i}\mathbf{x}\cdot\mathbf{x}} \varphi\left(-(\mathbf{u}-\mathbf{i}\mathbf{a})\right) \int\_{-\infty}^{\infty} e^{\mathrm{i}\mathbf{u}\cdot\mathbf{y}+\mathbf{a}\cdot\mathbf{y}} \upsilon(\mathbf{y},t\_{k+1})d\mathbf{y}d\mathbf{u}.\tag{26}$$

Denote

$$\kappa = \delta \sqrt{-\frac{\partial^2 \varrho\_T(\boldsymbol{u})}{\partial \boldsymbol{u}^2} \Big|\_{\boldsymbol{u}=\boldsymbol{0}} + \left(\frac{\partial \varrho\_T(\boldsymbol{u})}{\partial \boldsymbol{u}} \Big|\_{\boldsymbol{u}=\boldsymbol{0}}\right)^2}$$

where *ϕT*(*u*) is the characteristic function of *XT*, and *δ* is a proportionality constant. According to the suggestion from Lord *et al* (2008), we can take *δ* = 20 for the GBM model and *δ* = 40 for other exponential Lévy models. Let *N* be a power of 2. We consider the grid points on *x*-axes:

$$\mathbf{x}\_{j} = \left(j - \frac{1}{2}N\right)\Delta x, \quad j = 0, 1, \dots, N - 1, \dots$$

where Δ*x* = *κ*/*N*. Furthermore, we also consider the grid points for the numerical integrals in (26):

$$
\mu\_j = (j - \frac{1}{2}N)\Delta\mu\_\prime \quad y\_j = x\_{j\prime} \quad j = 0, 1, \ldots, N - 1, \ldots
$$

Here *S* is the spot price of underlying asset, *r* > 0 is the interest rate. Let *Xt* = log(*St*/*K*) be the logarithm of the underlying asset price *St* over the strike price *K*, and denote *x* = log(*S*/*K*) and *h* = log(*H*/*K*). Let *f*(· | *x*) be the condition density of *Xtk*<sup>+</sup><sup>1</sup> given *Xtk* = *x* for *tk* ∈ **T**. Set

Fourier Transform Methods for Option Pricing 283

*<sup>g</sup>*(*x*) = � *<sup>K</sup>*(*e<sup>x</sup>* <sup>−</sup> <sup>1</sup>)+, for a call option,

<sup>−</sup><sup>∞</sup> *<sup>v</sup>*(*y*, *tk*<sup>+</sup>1)*f*(*<sup>y</sup>* <sup>|</sup> *<sup>x</sup>*)*dy*, *tk* <sup>∈</sup> **<sup>T</sup>**,

Then the backward induction (23) and the price formula (30) can be rewritten by

*<sup>v</sup>*(*x*, *tk*) = *<sup>c</sup>*(*x*, *tk*))1{*x*<*h*} <sup>+</sup> *<sup>R</sup>*01{*x*≥*h*}, *tk* <sup>∈</sup> **<sup>T</sup>** \ **<sup>T</sup>***e*,

�

−*r*Δ*t* � ∞ −∞

We consider the infinite integrals *c*(*x*, *tk*) in (31). Since *f*(*y*|*x*) decays to zero very quickly as *y* → ±∞ we may choose two bounds *a* and *b*, which can be selected by (21), such that

without losing some significant accuracy. Note that the density *f*(*y* | *x*) has the following

<sup>−</sup>*r*Δ*<sup>t</sup>* <sup>∞</sup> ∑ *j*=0

*<sup>v</sup>*(*y*, *tk*<sup>+</sup>1) cos �

*<sup>f</sup>*(*<sup>y</sup>* <sup>|</sup> *<sup>x</sup>*) cos �

Since the Fourier-cosine expansion has a high accuracy with a few terms, we can truncate the

*N*−1 ∑ *j*=0

� � *<sup>b</sup> a*

<sup>2</sup> and *wj* = 1 for all *j* = 1, 2, 3, . . .. We substitute this expansion into the integral

*<sup>j</sup><sup>π</sup> <sup>y</sup>* <sup>−</sup> *<sup>a</sup> b* − *a*

<sup>2</sup> *<sup>e</sup>*

� *b a*

� *b a*

−*r*Δ*t*

*g*(*x*), *c*(*x*, *tk*)

*v*(*x*, *t*0) = *c*(*x*, *t*0) = *e*

−*r*Δ*t* � *b a*

> ∞ ∑ *j*=0 � *wj* cos �

*<sup>c</sup>*¯(*x*, *tk*) = (*<sup>b</sup>* <sup>−</sup> *<sup>a</sup>*)

*b* − *a*

*b* − *a*

<sup>2</sup> *<sup>e</sup>*

infinity series (34) and approximate *c*¯(*x*, *tk*) by leaving the first *N* terms, i.e.,

*Vj*(*tk*<sup>+</sup>1) = <sup>2</sup>

*Fj*(*x*) = <sup>2</sup>

*<sup>c</sup>*ˆ(*x*, *tk*) = (*<sup>b</sup>* <sup>−</sup> *<sup>a</sup>*)

*<sup>v</sup>*(*x*, *tML*) = *<sup>g</sup>*(*x*)1{*x*<*h*} <sup>+</sup> *<sup>R</sup>*01{*x*≥*h*},

*<sup>c</sup>*(*x*, *tk*) = *<sup>e</sup>*−*r*Δ*<sup>t</sup>* � <sup>∞</sup>

*<sup>v</sup>*(*x*, *tk*) = max �

where *<sup>v</sup>*(*x*, *tk*) = *<sup>V</sup>*(*Kex*, *tk*) for any *tk* <sup>∈</sup> **<sup>T</sup>**.

*c*¯(*x*, *tk*) = *e*

*b* − *a*

Fourier-cosine expansion on [*a*, *b*]:

(33) and then we can rewrite it by

where for each *j* = 1, 2, 3, . . .,

where *w*<sup>0</sup> = <sup>1</sup>

and

*<sup>f</sup>*(*<sup>y</sup>* <sup>|</sup> *<sup>x</sup>*) = <sup>2</sup>

⎧ ⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎩

and

*<sup>K</sup>*(<sup>1</sup> <sup>−</sup> *<sup>e</sup>x*)+, for a put option.

<sup>1</sup>{*x*<*h*} <sup>+</sup> *<sup>e</sup>*−*r*(*T*−*tk* )*R*01{*x*≥*h*}, *tk* <sup>∈</sup> **<sup>T</sup>***e*,

*v*(*y*, *tm*+1)*f*(*y* | *x*)*dy* ≈ *c*(*x*, *tk*), *tk* ∈ **T**. (33)

*jπ u* − *a b* − *a*

*wjFj*(*x*)*Vj*(*tk*<sup>+</sup>1), (34)

� *du*� ,

*dy*, (35)

*dy*. (36)

*<sup>f</sup>*(*<sup>u</sup>* <sup>|</sup> *<sup>x</sup>*) cos �

*<sup>j</sup><sup>π</sup> <sup>y</sup>* <sup>−</sup> *<sup>a</sup> b* − *a*

�

*<sup>j</sup><sup>π</sup> <sup>y</sup>* <sup>−</sup> *<sup>a</sup> b* − *a* �

*wjFj*(*x*)*Vj*(*tk*<sup>+</sup>1) ≈ *c*¯(*x*, *tk*), (37)

*v*(*y*, *t*1)*f*(*y* | *x*)*dy*, (32)

(31)

where Δ*u* = 2*π*/*κ*. It is clear that these grids satisfy the Nyquist relation: Δ*u*Δ*y* = 2*π*/*N*. Now, for each *tk* ∈ **T** and each *p* = 0, 1, . . . , *N*−1, approximating the integral in (26) with composite trapezoidal rule and the second integral with left rectangle rule yields

$$\mathcal{L}\left(\mathbf{x}\_{p},t\_{k}\right) \approx e^{-\alpha \mathbf{x}\_{p} - r\Delta t} \frac{\Delta \mathbf{u} \Delta \mathbf{y}}{2\pi} \sum\_{j=0}^{N-1} \left( e^{-\mathbf{i}\mathbf{u}\_{j}\mathbf{x}\_{p}} \boldsymbol{\varrho}\left(-\left(\mathbf{u}\_{j} - \mathbf{i}\mathbf{a}\right)\right) \sum\_{n=0}^{N-1} \omega\_{n} e^{\mathbf{i}\mathbf{u}\_{j}\mathbf{y}\_{n} + \mathbf{a}\mathbf{y}\_{n}} \boldsymbol{\upsilon}\left(\mathbf{y}\_{n},t\_{k+1}\right) \right)$$

where the weights *<sup>ω</sup><sup>n</sup>* are chosen as *<sup>ω</sup>*<sup>0</sup> <sup>=</sup> *<sup>ω</sup>N*−<sup>1</sup> <sup>=</sup> <sup>1</sup> <sup>2</sup> and *ω<sup>n</sup>* = 1 for *n* = 1, . . . , *N*−2. Noting that *<sup>u</sup>*<sup>0</sup> <sup>=</sup> <sup>−</sup><sup>1</sup> <sup>2</sup> *N*Δ*u* and Δ*x* = Δ*y* = 2*π*/(*N*Δ*u*), we have

$$c(\mathbf{x}\_{p^\*}t\_k) \approx e^{-a\mathbf{x}\_p - r\Delta t}e^{i\mathbf{l}\mathbf{u}\_0(y\_0 - \mathbf{x}\_0)}(-1)^p \sum\_{j=0}^{N-1} \left( e^{-i j p \frac{2\pi}{N}} e^{j j (y\_0 - \mathbf{x}\_0)\Delta t} \varrho\left(-(u\_j - \mathbf{i}\mathbf{a})\right) \right)$$

$$\cdot \frac{1}{N} \sum\_{n=0}^{N-1} e^{i j n \frac{2\pi}{N}} (-1)^n \omega\_n e^{\mathbf{a} y\_n} \upsilon(y\_n, t\_{k+1}) \, \Big|\,. \tag{27}$$

Now, we can employ the FFT algorithm to calculate the summations in the right side of (27). Once the integral *c*(*xp*, *tk*) is computed, we can determine the early-exercise price *S*<sup>∗</sup> *tk* , *tk* ∈ **T***e*, by the procedure: for each *tk* ∈ **T***e*, locate *jk* such that

$$\mathcal{L}(\mathbf{x}\_{\mathbf{j}\_{\ell}}, \mathbf{t}\_{k}) - \mathbf{g}(\mathbf{x}\_{\mathbf{j}\_{\ell}}) = \mathbf{0},\tag{28}$$

or,

$$\left(\mathcal{c}(\mathbf{x}\_{j\_k}, t\_k) - \mathbf{g}(\mathbf{x}\_{j\_k})\right) \left(\mathcal{c}(\mathbf{x}\_{j\_k+1}, t\_k) - \mathbf{g}(\mathbf{x}\_{j\_k+1})\right) \le 0. \tag{29}$$

In the case (28) set *<sup>x</sup>*∗(*tk*) = *xjk* , and in the case (29) set *<sup>x</sup>*∗(*tk*) = <sup>1</sup> <sup>2</sup> (*xjk* + *xjk*<sup>+</sup><sup>1</sup> ). Then the early-exercise price at every *tk* ∈ **T***<sup>e</sup>* is given by *S*<sup>∗</sup> *tk* <sup>=</sup> *<sup>S</sup>*0*ex*∗(*tk* ). Ding *et al* (2011b) gave a detail algorithm, which summarizes the above procedure, for pricing an up-and-out Bermudan barrier option, and the corresponding numerical experiments.

#### **4.2 The COS method for pricing Bermudan barrier options**

Recently, Fang & Oosterlee (2009) extended their COS method to price discrete early-exercise options under exponential Lévy models, and Fang & Oosterlee (2011) further considered such pricing problems under Heston's model.

Assume that the price process of the underlying asset *St* follows an exponential Lévy model, and **T** and **T***<sup>e</sup>* are the set of pre-specified monitored dates and the set of pre-specified exercise dates, respectively, before the maturity *T*, which are defined by (22). In the following, we apply the COS method to price the Bermudan barrier option which defined in preceding subsection, whose payoff is given by

$$\{G(S\_{l\_k})1\_{\{S\_{l\_k}$$

where *H* > *K* is the constant barrier and *R*<sup>0</sup> is the contractual rebate. Denote *V*(*S*, *tk*), *tk* ∈ **T***e*, the value of this Bermudan barrier option at time *tk* and the spot price *Stk* = *S*. As in the preceding subsection, with help of the risk-neutral valuation formula, this price process can be computed recursively by the backward induction (23). In specialty, the initial price is given by

$$V(\mathcal{S}, t\_0) = \mathcal{C}(\mathcal{S}, t\_0) = \mathbb{E}\left[e^{-r\Delta t} V(\mathcal{S}, t\_1) \mid \mathcal{S}\_{t\_0} = \mathcal{S}\right].\tag{30}$$

Here *S* is the spot price of underlying asset, *r* > 0 is the interest rate. Let *Xt* = log(*St*/*K*) be the logarithm of the underlying asset price *St* over the strike price *K*, and denote *x* = log(*S*/*K*) and *h* = log(*H*/*K*). Let *f*(· | *x*) be the condition density of *Xtk*<sup>+</sup><sup>1</sup> given *Xtk* = *x* for *tk* ∈ **T**. Set

$$g(\mathbf{x}) = \begin{cases} K(e^{\mathbf{x}} - 1)^{+}, & \text{for a call option,} \\ K(1 - e^{\mathbf{x}})^{+}, & \text{for a put option.} \end{cases}$$

Then the backward induction (23) and the price formula (30) can be rewritten by

$$\begin{cases} v(\mathbf{x}, t\_{ML}) = g(\mathbf{x}) \mathbf{1}\_{\{\mathbf{x} < \hbar\}} + R\_0 \mathbf{1}\_{\{\mathbf{x} \ge \hbar\}}, \\ c(\mathbf{x}, t\_k) \\ \begin{cases} c(\mathbf{x}, t\_k) \\ v(\mathbf{x}, t\_k) \end{cases} = e(\mathbf{x}, t\_k) \mathbf{1}\_{\{\mathbf{x} < \hbar\}} + R\_0 \mathbf{1}\_{\{\mathbf{x} \ge \hbar\}}, \quad t\_k \in \mathbb{T} \backslash \mathbb{T}\_{\mathcal{E}\_\mathcal{I}} \\ v(\mathbf{x}, t\_k) \quad = \max\left\{ g(\mathbf{x}), c(\mathbf{x}, t\_k) \right\} \mathbf{1}\_{\{\mathbf{x} < \hbar\}} + e^{-r(T - t\_k)} R\_0 \mathbf{1}\_{\{\mathbf{x} \ge \hbar\}^\mathcal{I}} \quad t\_k \in \mathbb{T}\_{\mathcal{E}\_\mathcal{I}} \end{cases} \tag{31}$$

and

18 Will-be-set-by-IN-TECH

where Δ*u* = 2*π*/*κ*. It is clear that these grids satisfy the Nyquist relation: Δ*u*Δ*y* = 2*π*/*N*. Now, for each *tk* ∈ **T** and each *p* = 0, 1, . . . , *N*−1, approximating the integral in (26) with

−(*uj*−i*α*)

 *N*−1 ∑ *n*=0

*ωne*

<sup>i</sup>*j*(*y*0−*x*0)Δ*uϕ*

*c*(*xjk* , *tk*) − *g*(*xjk* ) = 0, (28)

*tk* <sup>=</sup> *<sup>S</sup>*0*ex*∗(*tk* ). Ding *et al* (2011b) gave a detail

<sup>i</sup>*ujyn*+*αyn <sup>v</sup>*(*yn*, *tk*<sup>+</sup>1)

−(*uj*−i*α*)

. (27)

*tk*

≤ 0. (29)

<sup>2</sup> (*xjk* + *xjk*<sup>+</sup><sup>1</sup> ). Then the

, *tk* ∈ **T***e*,

<sup>2</sup> and *ω<sup>n</sup>* = 1 for *n* = 1, . . . , *N*−2. Noting

composite trapezoidal rule and the second integral with left rectangle rule yields

*N*−1 ∑ *j*=0

<sup>2</sup> *N*Δ*u* and Δ*x* = Δ*y* = 2*π*/(*N*Δ*u*), we have

i*u*0(*y*0−*x*0)

*<sup>N</sup>* (−1)*nωne*

(−1)*<sup>p</sup>*

Once the integral *c*(*xp*, *tk*) is computed, we can determine the early-exercise price *S*<sup>∗</sup>

*N*−1 ∑ *j*=0

Now, we can employ the FFT algorithm to calculate the summations in the right side of (27).

algorithm, which summarizes the above procedure, for pricing an up-and-out Bermudan

Recently, Fang & Oosterlee (2009) extended their COS method to price discrete early-exercise options under exponential Lévy models, and Fang & Oosterlee (2011) further considered such

Assume that the price process of the underlying asset *St* follows an exponential Lévy model, and **T** and **T***<sup>e</sup>* are the set of pre-specified monitored dates and the set of pre-specified exercise dates, respectively, before the maturity *T*, which are defined by (22). In the following, we apply the COS method to price the Bermudan barrier option which defined in preceding

*<sup>k</sup>* <sup>&</sup>lt;*H*} <sup>+</sup> *<sup>R</sup>*01{*St*

where *H* > *K* is the constant barrier and *R*<sup>0</sup> is the contractual rebate. Denote *V*(*S*, *tk*), *tk* ∈ **T***e*, the value of this Bermudan barrier option at time *tk* and the spot price *Stk* = *S*. As in the preceding subsection, with help of the risk-neutral valuation formula, this price process can be computed recursively by the backward induction (23). In specialty, the initial price is given

> *e* −*r*Δ*t*

*<sup>k</sup>* <sup>≥</sup>*H*}, *tk* <sup>∈</sup> **<sup>T</sup>**,

*V*(*S*, *t*1) | *St*<sup>0</sup> = *S*

. (30)

 *e* <sup>−</sup>i*jp* <sup>2</sup>*<sup>π</sup> <sup>N</sup> e*

*<sup>α</sup>yn <sup>v</sup>*(*yn*, *tk*<sup>+</sup>1)

*<sup>c</sup>*(*xjk*+1, *tk*) <sup>−</sup> *<sup>g</sup>*(*xjk*<sup>+</sup>1)

*e*

*c*(*xjk* , *tk*) − *g*(*xjk* )

barrier option, and the corresponding numerical experiments.

*<sup>G</sup>*(*Stk* )1{*St*

*V*(*S*, *t*0) = *C*(*S*, *t*0) = **E**

**4.2 The COS method for pricing Bermudan barrier options**

In the case (28) set *<sup>x</sup>*∗(*tk*) = *xjk* , and in the case (29) set *<sup>x</sup>*∗(*tk*) = <sup>1</sup>

 *e* <sup>−</sup>i*ujxp ϕ* 

*c*(*xp*, *tk*) ≈ *e*

*c*(*xp*, *tk*) ≈ *e*

that *<sup>u</sup>*<sup>0</sup> <sup>=</sup> <sup>−</sup><sup>1</sup>

or,

by

<sup>−</sup>*αxp*−*r*Δ*<sup>t</sup>* <sup>Δ</sup>*u*Δ*<sup>y</sup>*

where the weights *<sup>ω</sup><sup>n</sup>* are chosen as *<sup>ω</sup>*<sup>0</sup> <sup>=</sup> *<sup>ω</sup>N*−<sup>1</sup> <sup>=</sup> <sup>1</sup>

−*αxp*−*r*Δ*t*

*N*−1 ∑ *n*=0 *e* i*jn* <sup>2</sup>*<sup>π</sup>*

by the procedure: for each *tk* ∈ **T***e*, locate *jk* such that

early-exercise price at every *tk* ∈ **T***<sup>e</sup>* is given by *S*<sup>∗</sup>

pricing problems under Heston's model.

subsection, whose payoff is given by

· 1 *N* 2*π*

$$v(\mathbf{x}, t\_0) = c(\mathbf{x}, t\_0) = e^{-r\Delta t} \int\_{-\infty}^{\infty} v(y, t\_1) f(y \mid \mathbf{x}) dy,\tag{32}$$

where *<sup>v</sup>*(*x*, *tk*) = *<sup>V</sup>*(*Kex*, *tk*) for any *tk* <sup>∈</sup> **<sup>T</sup>**.

We consider the infinite integrals *c*(*x*, *tk*) in (31). Since *f*(*y*|*x*) decays to zero very quickly as *y* → ±∞ we may choose two bounds *a* and *b*, which can be selected by (21), such that

$$\varepsilon(\mathbf{x}, t\_k) = e^{-r\Delta t} \int\_a^b v(y, t\_{m+1}) f(y \mid \mathbf{x}) dy \approx c(\mathbf{x}, t\_k), \quad t\_k \in \mathbb{T}.\tag{33}$$

without losing some significant accuracy. Note that the density *f*(*y* | *x*) has the following Fourier-cosine expansion on [*a*, *b*]:

$$f(y \mid \mathbf{x}) = \frac{2}{b-a} \sum\_{j=0}^{\infty} \left[ w\_j \cos \left( j \pi \frac{y-a}{b-a} \right) \int\_a^b f(u \mid \mathbf{x}) \cos \left( j \pi \frac{u-a}{b-a} \right) du \right] \mathbf{x}$$

where *w*<sup>0</sup> = <sup>1</sup> <sup>2</sup> and *wj* = 1 for all *j* = 1, 2, 3, . . .. We substitute this expansion into the integral (33) and then we can rewrite it by

$$\bar{c}(\mathbf{x}, t\_k) = \frac{(b - a)}{2} e^{-r\Delta t} \sum\_{j=0}^{\infty} w\_j F\_j(\mathbf{x}) V\_j(t\_{k+1}) \,\tag{34}$$

where for each *j* = 1, 2, 3, . . .,

$$V\_j(t\_{k+1}) = \frac{2}{b-a} \int\_a^b v(y, t\_{k+1}) \cos \left( j \pi \frac{y-a}{b-a} \right) dy,\tag{35}$$

and

$$F\_{\hat{f}}(\mathbf{x}) = \frac{2}{b-a} \int\_{a}^{b} f(y \mid \mathbf{x}) \cos \left( j\pi \frac{y-a}{b-a} \right) dy. \tag{36}$$

Since the Fourier-cosine expansion has a high accuracy with a few terms, we can truncate the infinity series (34) and approximate *c*¯(*x*, *tk*) by leaving the first *N* terms, i.e.,

$$\mathfrak{E}(\mathbf{x}, t\_k) = \frac{(b-a)}{2} e^{-r\Delta t} \sum\_{j=0}^{N-1} w\_j F\_j(\mathbf{x}) V\_j(t\_{k+1}) \approx \overline{\mathfrak{c}}(\mathbf{x}, t\_k), \tag{37}$$

We first calculate *Vj*(*tML*), *j* = 0, 1, . . . , *N*−1. We have

*<sup>e</sup><sup>x</sup>* cos

*jπ x* − *a b* − *a*

1 + ( *<sup>j</sup><sup>π</sup> <sup>b</sup>*−*<sup>a</sup>* )<sup>2</sup>

<sup>+</sup> *<sup>j</sup><sup>π</sup> b* − *a*

<sup>Ψ</sup>*j*(*x*1, *<sup>x</sup>*2) =

 cos *jπ x*<sup>2</sup> − *a b* − *a*

sin *jπ x*<sup>2</sup> − *a b* − *a*

sin *jπ x*<sup>2</sup> − *a b* − *a*

*Dj*(*x*1, *<sup>x</sup>*2) = <sup>2</sup>*R*<sup>0</sup>

*Cj*(*x*1, *<sup>x</sup>*2; *tML*) = <sup>2</sup>

*Vj*(*tk*) = *Cj*(*a*, *h*; *tk*) + *e*

for all *j* = 0, 1, . . . , *N*−1. Denote

*x*1

<sup>Φ</sup>*j*(*x*1, *<sup>x</sup>*2) = <sup>1</sup>

0, . . . , *N*−1, we have the following result:

admit the following analytic solutions:

<sup>Φ</sup>*j*(*x*1, *<sup>x</sup>*2) = *<sup>x</sup>*<sup>2</sup>

and

point *x*∗

i.e., *c*˜(*x*∗

*<sup>k</sup>* , *tk*) = *g*(*x*<sup>∗</sup>

secant method to find the root *x*∗

cases for an up-and-out barrier option:

*<sup>k</sup>* ). Let

Then, the problem becomes to find the root *x*∗

*Vj*(*tML*) = *Cj*(0, *<sup>h</sup>*; *tML*) + *Dj*(*h*, *<sup>b</sup>*), for a call option,

Fourier Transform Methods for Option Pricing 285

for any *a* ≤ *x*<sup>1</sup> ≤ *x*<sup>2</sup> ≤ *b* and *j* = 0, . . . , *N*−1. Then, by simple integration, these integrals

for *j* = 0, . . . , *N*−1, with Ψ0(*x*1, *x*2) = *x*<sup>2</sup> − *x*1. Moreover, by a simple calculation, *j* =

*b* − *a*

*b* − *a αK* 

where *α* is a parameter such that *α* = 1 for a call option and *α* = −1 for a put option.

−*r*(*T*−*tk*−<sup>1</sup>)

Since the integral *Dj*(*x*1, *x*2) has the analytic representation (45), we only need to calculate the integral *Cj*(*x*1, *x*2; *tk*). Fang & Oosterlee (2009) gave an efficient numerical algorithm which approximates *Cj*(*x*1, *x*2; *tk*) by using FFT method with the operation cost *O*(*N* log2(*N*)).

Finally, we consider to calculate the integrals *Vk*(*tk*) for *tk* ∈ **T***e*. It is clear that we should find the value *v*˜(*x*, *tk*) in the last equation in (31), or equivalently, to determine the early-exercise

*hk*(*y*) = *c*˜(*y*, *tk*) − *g*(*y*), *tk* ∈ **T***e*.

function *c*˜(*y*, *tk*), which is given in (40), is bounded and smooth, and the function *g*(*y*) is smooth except for *y* = 0 and bounded in [*a*, *b*]. We can use the Newton's method or the

*<sup>k</sup>* . Here if *x*<sup>∗</sup>

boundary point *a* or *b*. Once we find the early-exercise point *x*∗

*<sup>k</sup>* at each time *tk*, which is the point where the continuation value is equal to the payoff,

And next, we consider to calculate the integrals *Vj*(*tk*) for *tk* ∈ **T** \ **T***e*. We have

*Cj*(*a*, 0; *tML*) + *Dj*(*h*, *b*), for a put option,

*dx* and <sup>Ψ</sup>*k*(*x*1, *<sup>x</sup>*2) = *<sup>x</sup>*<sup>2</sup>

*<sup>e</sup>x*<sup>2</sup> <sup>−</sup> *<sup>j</sup><sup>π</sup> b* − *a*

<sup>−</sup> sin *jπ x*<sup>1</sup> − *a b* − *a*

*<sup>e</sup>x*<sup>2</sup> <sup>−</sup> cos

Φ*j*(*x*1, *x*2) − Ψ*j*(*x*1, *x*2)

*Dj*(*h*, *b*), *j* = 0, 1, . . . , *N*−1.

*x*1

*jπ x*<sup>1</sup> − *a b* − *a*

 *b* − *a <sup>j</sup><sup>π</sup>* ,

Ψ*j*(*x*1, *x*2), (42)

*<sup>k</sup>* of each equation *hk*(*y*) = 0. Note that the

*<sup>k</sup>* , *tk* ∈ **T***e*, we have two different

*<sup>k</sup>* is not in the interval [*a*, *b*], we set it in the nearest

sin *jπ x*<sup>1</sup> − *a b* − *a*

cos *jπ x* − *a b* − *a*

> *ex*1

> > *ex*1 ,

 *dx*,

, (43)

for any *tk* <sup>∈</sup> **<sup>T</sup>**, where *<sup>w</sup>*<sup>0</sup> <sup>=</sup> *wN*−<sup>1</sup> <sup>=</sup> <sup>1</sup> <sup>2</sup> and *wj* = 1 for *j* = 1, . . . , *N*−2. On the other hand, we can represent each integral (36) by approximation as the following

$$\begin{split} F\_{\hat{l}}(\mathbf{x}) &= \frac{2}{b-a} \text{Re} \Big\{ \exp \left( -\mathbf{i} \frac{ja\pi}{b-a} \right) \int\_{a}^{b} f(u \mid \mathbf{x}) \exp \left( \mathbf{i} \frac{ju\pi}{b-a} \right) du \Big\} \\ &\approx \frac{2}{b-a} \text{Re} \Big\{ \exp \left( -\mathbf{i} \frac{ja\pi}{b-a} \right) \int\_{-\infty}^{\infty} f(u \mid \mathbf{x}) \exp \left( \mathbf{i} \frac{j\pi}{b-a} u \right) du \Big\}. \end{split}$$

Let *ϕ*(*u*; *x*) be the characteristic function of *f*(· | *x*). Then, we can approximate each *Fj*(*x*) by

$$F\_{\vec{\}}(\mathbf{x}) \approx \frac{2}{b-a} \text{Re}\{ \exp\left( -\mathbf{i} \frac{ja\pi}{b-a} \right) \cdot \boldsymbol{\varrho} \left( \frac{j\pi}{b-a}; \mathbf{x} \right) \}.$$

And so, we get the further numerical approximations of the integrals *c*ˆ(*x*, *tk*) in (37) by

$$\mathfrak{E}(\mathbf{x}, t\_k) = e^{-r\Delta t} \sum\_{j=0}^{N-1} w\_j \text{Re}\left\{ \exp\left(-\mathbf{i}\frac{ja\pi}{b-a}\right) \cdot \boldsymbol{\varrho}\left(\frac{j\pi}{b-a}; \mathbf{x}\right) \right\} V\_j(t\_{k+1}) \approx \mathfrak{E}(\mathbf{x}, t\_j). \tag{38}$$

In special, the approximation of initial price *v*(*x*, *t*0) in (32) is given by

$$\vec{w}(\mathbf{x}, t\_0) = e^{-r\Delta t} \sum\_{j=0}^{N-1} w\_j \text{Re}\left\{ \exp\left(-\mathbf{i}\frac{ja\pi}{b-a}\right) \cdot \boldsymbol{\varrho}\left(\frac{j\pi}{b-a}; \mathbf{x}\right) \right\} V\_j(t\_1). \tag{39}$$

Meanwhile, since the characteristic function *φ*(*u*; *x*) possesses the property:

$$\varphi(u;\mathfrak{x}) = \mathfrak{q}(\mathfrak{u}) \cdot e^{\mathrm{i}u\mathfrak{x}}, \quad \mathfrak{u} \in \mathbb{R}\_{\mathsf{A}}$$

where *ϕ*(*u*) = *ϕ*(*u*; 0), the approximations (38) and (39) can be simplified to

$$\mathfrak{z}(\mathbf{x}, t\_k) = e^{-r\Lambda t} \sum\_{j=0}^{N-1} w\_j \text{Re}\left\{ \exp\left( \mathbf{i} j \pi \frac{\mathbf{x} - a}{b - a} \right) \cdot \boldsymbol{\varrho} \left( \frac{j \pi}{b - a} \right) \right\} V\_j(t\_{k+1}), \tag{40}$$

and the initial price of option

$$\vec{v}(\mathbf{x}, t\_0) = e^{-r\Delta t} \sum\_{j=0}^{N-1} w\_j \text{Re}\left\{ \exp\left(\mathbf{i}j\pi \frac{\mathbf{x} - a}{b - a}\right) \cdot \boldsymbol{\varrho}\left(\frac{j\pi}{b - a}\right) \right\} V\_j(t\_1). \tag{41}$$

In order to use this approximate formulation we still need to compute the integrals *Vj*(*tk*). For convenience we introduce the following notions: for any *a* ≤ *x*<sup>1</sup> ≤ *x*<sup>2</sup> ≤ *b*,

$$\mathbf{C}\_{j}(\mathbf{x}\_{1}, \mathbf{x}\_{2}; t\_{k}) = \frac{2}{b-a} \int\_{\mathcal{X}\_{1}}^{\mathbf{x}\_{2}} \tilde{c}(\mathbf{x}, t\_{k}) \cos \left( j\pi \frac{\mathbf{x} - a}{b - a} \right) d\mathbf{x}, \quad j = 0, 1, \ldots, N - 1, \ldots$$

and

$$D\_j(\mathbf{x}\_1, \mathbf{x}\_2) = \frac{2R\_0}{b-a} \int\_{\mathbf{x}\_1}^{\mathbf{x}\_2} \cos \left( j\pi \frac{\mathbf{x}-a}{b-a} \right) d\mathbf{x}\_1 \quad j = 0, 1, \dots, N-1, \quad 1$$

where *c*˜(*x*, *tML*) = *g*(*x*), and *c*˜(*x*, *tk*), *tk* ∈ **T**, are given in (40).

We first calculate *Vj*(*tML*), *j* = 0, 1, . . . , *N*−1. We have

$$V\_j(t\_{ML}) = \begin{cases} \mathbb{C}\_j(0, h; t\_{ML}) + D\_j(h, b), & \text{for a call option,} \\ \mathbb{C}\_j(a, 0; t\_{ML}) + D\_j(h, b), & \text{for a put option,} \end{cases}$$

for all *j* = 0, 1, . . . , *N*−1. Denote

$$\Phi\_j(\mathbf{x}\_1, \mathbf{x}\_2) = \int\_{\mathbf{x}\_1}^{\mathbf{x}\_2} e^{\mathbf{x}} \cos \left( j \pi \frac{\mathbf{x} - a}{b - a} \right) d\mathbf{x} \quad \text{and} \quad \Psi\_k(\mathbf{x}\_1, \mathbf{x}\_2) = \int\_{\mathbf{x}\_1}^{\mathbf{x}\_2} \cos \left( j \pi \frac{\mathbf{x} - a}{b - a} \right) d\mathbf{x} \,\mathrm{d}\varphi$$

for any *a* ≤ *x*<sup>1</sup> ≤ *x*<sup>2</sup> ≤ *b* and *j* = 0, . . . , *N*−1. Then, by simple integration, these integrals admit the following analytic solutions:

$$\begin{split} \Phi\_{j}(\mathbf{x}\_{1}, \mathbf{x}\_{2}) &= \frac{1}{1 + (\frac{j\pi}{b-a})^{2}} \Big[ \cos \left( j\pi \frac{\mathbf{x}\_{2} - a}{b - a} \right) e^{\mathbf{x}\_{2}} - \cos \left( j\pi \frac{\mathbf{x}\_{1} - a}{b - a} \right) e^{\mathbf{x}\_{1}} \\ &+ \frac{j\pi}{b - a} \sin \left( j\pi \frac{\mathbf{x}\_{2} - a}{b - a} \right) e^{\mathbf{x}\_{2}} - \frac{j\pi}{b - a} \sin \left( j\pi \frac{\mathbf{x}\_{1} - a}{b - a} \right) e^{\mathbf{x}\_{1}} \Big]. \end{split}$$

and

20 Will-be-set-by-IN-TECH

 *<sup>b</sup> a*

 <sup>∞</sup> −∞

<sup>−</sup> <sup>i</sup> *ja<sup>π</sup> b* − *a*

> · *ϕ jπ b* − *a* ; *x*

<sup>−</sup> <sup>i</sup> *ja<sup>π</sup> b* − *a*  · *ϕ jπ b* − *a* ; *x* 

<sup>i</sup>*ux*, *<sup>u</sup>* <sup>∈</sup> **<sup>R</sup>**,

 · *ϕ jπ b* − *a*

> · *ϕ jπ b* − *a*

*dx*, *j* = 0, 1, . . . , *N*−1,

*dx*, *j* = 0, 1, . . . , *N*−1,

 · *ϕ jπ b* − *a* ; *x* .

Let *ϕ*(*u*; *x*) be the characteristic function of *f*(· | *x*). Then, we can approximate each *Fj*(*x*) by

And so, we get the further numerical approximations of the integrals *c*ˆ(*x*, *tk*) in (37) by

<sup>−</sup> <sup>i</sup> *ja<sup>π</sup> b* − *a* *f*(*u* | *x*) exp

*f*(*u* | *x*) exp

<sup>−</sup> <sup>i</sup> *ja<sup>π</sup> b* − *a*

<sup>−</sup> <sup>i</sup> *ja<sup>π</sup> b* − *a*

<sup>2</sup> and *wj* = 1 for *j* = 1, . . . , *N*−2. On the other hand, we

 <sup>i</sup> *ju<sup>π</sup> b* − *a*

 <sup>i</sup> *<sup>j</sup><sup>π</sup> b* − *a u du* .

 *du* 

*Vj*(*tk*<sup>+</sup>1) ≈ *c*ˆ(*x*, *tj*). (38)

*Vj*(*t*1). (39)

*Vj*(*tk*<sup>+</sup>1), (40)

*Vj*(*t*1). (41)

for any *tk* <sup>∈</sup> **<sup>T</sup>**, where *<sup>w</sup>*<sup>0</sup> <sup>=</sup> *wN*−<sup>1</sup> <sup>=</sup> <sup>1</sup>

*Fj*(*x*) = <sup>2</sup>

*c*˜(*x*, *tk*) = *e*

*b* − *a* Re exp 

*Fj*(*x*) <sup>≈</sup> <sup>2</sup>

*N*−1 ∑ *j*=0

−*r*Δ*t*

−*r*Δ*t*

−*r*Δ*t*

*b* − *a*

*b* − *a*

where *c*˜(*x*, *tML*) = *g*(*x*), and *c*˜(*x*, *tk*), *tk* ∈ **T**, are given in (40).

*Dj*(*x*1, *<sup>x</sup>*2) = <sup>2</sup>*R*<sup>0</sup>

*N*−1 ∑ *j*=0

> *N*−1 ∑ *j*=0

*b* − *a* Re exp 

*wj*Re exp 

In special, the approximation of initial price *v*(*x*, *t*0) in (32) is given by

*wj*Re exp 

Meanwhile, since the characteristic function *φ*(*u*; *x*) possesses the property:

where *ϕ*(*u*) = *ϕ*(*u*; 0), the approximations (38) and (39) can be simplified to

*wj*Re exp i*jπ x* − *a b* − *a*

*wj*Re exp i*jπ x* − *a b* − *a*

convenience we introduce the following notions: for any *a* ≤ *x*<sup>1</sup> ≤ *x*<sup>2</sup> ≤ *b*,

 *x*<sup>2</sup> *x*1 cos *jπ x* − *a b* − *a*

 *x*<sup>2</sup> *x*1

In order to use this approximate formulation we still need to compute the integrals *Vj*(*tk*). For

 *jπ x* − *a b* − *a*

*c*˜(*x*, *tk*) cos

*ϕ*(*u*; *x*) = *ϕ*(*u*) · *e*

*N*−1 ∑ *j*=0

<sup>≈</sup> <sup>2</sup> *b* − *a* Re exp 

−*r*Δ*t*

*v*˜(*x*, *t*0) = *e*

*c*˜(*x*, *tk*) = *e*

*v*˜(*x*, *t*0) = *e*

*Cj*(*x*1, *<sup>x</sup>*2; *tk*) = <sup>2</sup>

and the initial price of option

and

can represent each integral (36) by approximation as the following

$$\Psi\_j(\mathbf{x}\_1, \mathbf{x}\_2) = \left[ \sin \left( j \pi \frac{\mathbf{x}\_2 - a}{b - a} \right) - \sin \left( j \pi \frac{\mathbf{x}\_1 - a}{b - a} \right) \right] \frac{b - a}{j \pi} \lambda$$

for *j* = 0, . . . , *N*−1, with Ψ0(*x*1, *x*2) = *x*<sup>2</sup> − *x*1. Moreover, by a simple calculation, *j* = 0, . . . , *N*−1, we have the following result:

$$D\_j(\mathbf{x}\_1, \mathbf{x}\_2) = \frac{2\mathbf{R}\_0}{b-a} \Psi\_j(\mathbf{x}\_1, \mathbf{x}\_2),\tag{42}$$

$$\mathcal{C}\_{j}(\mathbf{x}\_{1}, \mathbf{x}\_{2}; t\_{ML}) = \frac{2}{b-a} aK \{\Phi\_{j}(\mathbf{x}\_{1}, \mathbf{x}\_{2}) - \Psi\_{j}(\mathbf{x}\_{1}, \mathbf{x}\_{2})\},\tag{43}$$

where *α* is a parameter such that *α* = 1 for a call option and *α* = −1 for a put option.

And next, we consider to calculate the integrals *Vj*(*tk*) for *tk* ∈ **T** \ **T***e*. We have

$$V\_j(t\_k) = \mathbb{C}\_j(a, b; t\_k) + e^{-r(T - t\_{k-1})} D\_j(h, b), \quad j = 0, 1, \dots, N - 1.$$

Since the integral *Dj*(*x*1, *x*2) has the analytic representation (45), we only need to calculate the integral *Cj*(*x*1, *x*2; *tk*). Fang & Oosterlee (2009) gave an efficient numerical algorithm which approximates *Cj*(*x*1, *x*2; *tk*) by using FFT method with the operation cost *O*(*N* log2(*N*)).

Finally, we consider to calculate the integrals *Vk*(*tk*) for *tk* ∈ **T***e*. It is clear that we should find the value *v*˜(*x*, *tk*) in the last equation in (31), or equivalently, to determine the early-exercise point *x*∗ *<sup>k</sup>* at each time *tk*, which is the point where the continuation value is equal to the payoff, i.e., *c*˜(*x*∗ *<sup>k</sup>* , *tk*) = *g*(*x*<sup>∗</sup> *<sup>k</sup>* ). Let

$$h\_k(y) = \tilde{c}(y, t\_k) - g(y), \quad t\_k \in \mathbb{T}\_{\mathcal{C}}.$$

Then, the problem becomes to find the root *x*∗ *<sup>k</sup>* of each equation *hk*(*y*) = 0. Note that the function *c*˜(*y*, *tk*), which is given in (40), is bounded and smooth, and the function *g*(*y*) is smooth except for *y* = 0 and bounded in [*a*, *b*]. We can use the Newton's method or the secant method to find the root *x*∗ *<sup>k</sup>* . Here if *x*<sup>∗</sup> *<sup>k</sup>* is not in the interval [*a*, *b*], we set it in the nearest boundary point *a* or *b*. Once we find the early-exercise point *x*∗ *<sup>k</sup>* , *tk* ∈ **T***e*, we have two different cases for an up-and-out barrier option:

Thus, letting *Pt*<sup>1</sup> *v*(*x*) = **E**

induction (44) we get

where *vk*(*x*) = *er*Δ*<sup>t</sup>*

Note that 1{*x*>0} <sup>=</sup> <sup>1</sup>

transform) we obtain

and hence,

�

⎧ ⎪⎨

⎪⎩

1 + sgn(*x*)

from the backward induction (45).

<sup>1</sup>{*x*>*a*} · *<sup>v</sup>*(*x*) = <sup>1</sup>

and using equation (46) we obtain

*Xtk* = *x* possesses the property:

*Pt*<sup>1</sup> *vk*<sup>+</sup>1(*x*) in (45) becomes to

we have

F�

*Pt*<sup>1</sup> *vk*(*x*) �

F�

F�

2 �

<sup>1</sup>{*x*>*a*} · *<sup>v</sup>*

<sup>1</sup>{*x*>*a*} · *<sup>v</sup>*

*Pt*<sup>1</sup> *vk*<sup>+</sup>1(*x*) = � <sup>+</sup><sup>∞</sup>

(*u*) = <sup>F</sup>�

2 � *v*(*Xt*<sup>1</sup> ) | *Xt*<sup>0</sup> = *x*

*v*0(*x*) = *Pt*<sup>1</sup> *v*1(*x*),

�

<sup>F</sup>(1{*x*>0} · *<sup>v</sup>*)(*u*) = <sup>1</sup>

Let *a* ∈ **R** and T*<sup>a</sup>* be the transform operator: T*av*(*x*) = *v*(*x* − *a*). Then, we have

<sup>1</sup>{*x*>*a*} <sup>=</sup> <sup>T</sup>*a*1{*x*>0} <sup>=</sup> <sup>1</sup>

*v*(*x*) + *v*(*x*) · T*a*sgn(*x*)

Taking the Fourier transform to both sides of this equation, we have

� (*u*) = <sup>1</sup> 2

� (*u*) = <sup>1</sup> 2

−∞

*vM*(*x*) = <sup>1</sup>{*x*>*l*} · *<sup>K</sup>*(<sup>1</sup> <sup>−</sup> *<sup>e</sup>x*)+,

�

Fourier Transform Methods for Option Pricing 287

*vk*(*x*) = <sup>1</sup>{*x*>*l*}*Pt*<sup>1</sup> *vk*<sup>+</sup>1(*x*), *<sup>j</sup>* <sup>=</sup> *<sup>M</sup>*−1, . . . , 1,

2 �

> 2 �

> > � <sup>=</sup> <sup>1</sup> 2 �

2 F � T*a* �

2 i*e* <sup>i</sup>*au*H� *e*

> � +∞ −∞

> > *f*(−*x*) � (*u*) · F�

<sup>F</sup>*v*(*u*) + <sup>1</sup>

<sup>F</sup>*v*(*u*) + <sup>1</sup>

On other hand, noting that, for each *k* = 1, . . . , *M*−1, the condition density of *Xtk*<sup>+</sup><sup>1</sup> given

*f*(*y*|*x*) = *f*(*y*−*x*), *x*, *y* ∈ **R**, where *f*(*y*) is the density of *Xt*<sup>1</sup> under the initial condition *Xt*<sup>0</sup> = *x*. The infinite integrals

*vk*<sup>+</sup>1(*y*)*f*(*y*−*x*)d*y* =

*f*(−*x*) ∗ *vk*<sup>+</sup>1(*x*)

Then, this integral can be rewritten as a convolution of *vk*<sup>+</sup>1(*x*+*z*) and the function *f*(−*x*), i.e. *Pt*<sup>1</sup> *vk*<sup>+</sup>1(*x*) = *f*(−*x*) ∗ *vk*<sup>+</sup>1(*x*). Note that supp*vk*(*x*) ⊂ (*l*, 0) for each *k* = *M*, . . . , 1 from the backward induction (45). We can take the Fourier transform in the both sides of (45). Applying the Property P3 (Convolution)

�

(*u*) = <sup>F</sup>�

*V*(*x*, *tk*). And hence, the problem now becomes to find the function *v*0(*x*)

F*v*(*u*) + iH(F*v*)(*u*)

1 + T*a*sgn(*x*)

for all *x* ∈ **R**. Using the Property P4 (Relation with Hilbert

� ,

�

�� (*u*).

�

*vk*<sup>+</sup>1(*x*+*z*)*f*(*z*)*dz*.

*vk*<sup>+</sup>1(*x*) � (*u*).

sgn · T−*av*

� (*x*) � .

(*u*). (47)

*v*(*x*) + T*<sup>a</sup>*

sgn · T−*av*

<sup>−</sup>i*ay*F*v*(*y*)

�

. (46)

be the expectation operator, from the backward

(45)

Case 1: *x*∗ *<sup>k</sup>* < *h*, which means the early-exercise point doesn't hit the up barrier. Thus, We have the authority to decide to execute the option now or reserve it to the next time point. So we can split the integral that defines *Vj*(*tk*) into three parts: [*a*, *x*<sup>∗</sup> *<sup>k</sup>* ], (*x*<sup>∗</sup> *<sup>k</sup>* , *h*) and [*h*, *b*]. We have

$$V\_{\boldsymbol{j}}(t\_k) = \begin{cases} \mathbb{C}\_{\boldsymbol{j}}(\boldsymbol{a}, \mathbf{x}\_k^\*; t\_k) + \mathbb{C}\_{\boldsymbol{j}}(\mathbf{x}\_{k'}^\* \boldsymbol{h}; t\_{ML}) + e^{-r(T - t\_{k-1})} \boldsymbol{D}\_{\boldsymbol{j}}(\boldsymbol{h}, \boldsymbol{b}), \text{ for a call option,} \\\mathbb{C}\_{\boldsymbol{j}}(\boldsymbol{a}, \mathbf{x}\_{k'}^\*; t\_{ML}) + \mathbb{C}\_{\boldsymbol{j}}(\mathbf{x}\_{k'}^\* \boldsymbol{h}; t\_k) + e^{-r(T - t\_{k-1})} \boldsymbol{D}\_{\boldsymbol{j}}(\boldsymbol{h}, \boldsymbol{b}), \text{ for a put option.} \end{cases}$$

Case 2: *x*∗ *<sup>k</sup>* ≥ *h* , which means the early-exercise point hits the up barrier. Thus, the option integral can be split into two parts: [*a*, *h*) and [*h*, *b*]:

$$V\_{\boldsymbol{j}}(t\_k) = \begin{cases} \mathbb{C}\_{\boldsymbol{j}}(\boldsymbol{a}, \boldsymbol{h}; t\_k) + e^{-r(T - t\_{k-1})} D\_{\boldsymbol{j}}(\boldsymbol{h}, \boldsymbol{b}), & \text{for a call option,} \\\mathbb{C}\_{\boldsymbol{j}}(\boldsymbol{a}, \boldsymbol{h}; t\_{ML}) + e^{-r(T - t\_{k-1})} D\_{\boldsymbol{j}}(\boldsymbol{h}, \boldsymbol{b}), & \text{for a put option,} \end{cases}$$

Ding *et al* (2011a) gave a detail algorithm, which summarizes the above procedure, for pricing an up-and-out Bermudan barrier option, and the corresponding numerical experiments.

#### **4.3 The fast Hilbert transform approach for pricing barrier options**

Feng & Linetsky (2008) presented a new numerical method to price discretely monitored barrier options under exponential Lévy models. Their method involves the relation with Hilbert transform (Property P4) and the Sine expansion in Hardy spaces. They also gave an efficient computational algorithm based on the fast Hilbert transform.

Let **T** = {*tk* : *k* = 1, . . . , *M*} be the set of pre-specified monitored dates, where

$$0 = t\_0 < t\_1 < \dots < t\_M = T \quad \text{with} \quad \Delta t = t\_k - t\_{k-1} = T/M.$$

We consider a European-type barrier put option whose payoff at maturity *T* is given by

$$G = \mathbf{1}\_{\{S\_{l\_1} > L\}} \mathbf{1}\_{\{S\_{l\_2} > L\}} \cdots \mathbf{1}\_{\{S\_{l\_M} > L\}} (\mathbf{K} - \mathbf{S\_T})^+ \ . $$

where *K* is the strike price and 0 < *L* < *K* is the lower barrier. We also assume that the underlying asset price process is given by *St* = *KeXt* with *S*<sup>0</sup> = *S*, where *Xt* is a Lévy process started at *x* = log(*S*/*K*). Denoting *l* = log(*L*/*K*), then, with help of the risk-neutral valuation formula, the price of this option is given by

$$V(\mathbf{x}, t\_0) = e^{-rT} \mathbb{E} \left[ \mathbf{1}\_{\{X\_{l\_1} > l\}} \mathbf{1}\_{\{X\_{l\_2} > l\}} \cdot \dots \mathbf{1}\_{\{X\_{l\_M} > l\}} K (1 - e^{X\_{l\_M}})^+ \mid X\_0 = \mathbf{x} \right]\_{\mathcal{X}}$$

which can be computed recursively by the following backward induction:

$$\begin{cases} V(\mathbf{x}, t\_M) = \mathbf{1}\_{\{\mathbf{x} > l\}} K (1 - e^{\mathbf{x}})^+, \\ V(\mathbf{x}, t\_k) \quad = e^{-r\Lambda t} \mathbf{1}\_{\{\mathbf{x} > l\}} \mathbb{E} \left[ V(\mathbf{X}\_{t\_{l+1}}, t\_{k+1}) \mid \mathbf{X}\_{t\_k} = \mathbf{x} \right], \quad k = M - 1, \dots, 1, \\ V(\mathbf{x}, t\_0) \quad = e^{-rT} \mathbb{E} \left[ e^{-r\Lambda t} V(\mathbf{X}\_{t\_1}, t\_1) \mid \mathbf{X}\_{t\_0} = \mathbf{x} \right]. \end{cases} \tag{44}$$

Since each Lévy process is stationary and has independent increments, for each *k* = 1, . . . , *M*−1, we have

$$\mathbb{E}\left[V(\mathbf{X}\_{t\_{k+1}}, t\_{k+1}) \mid \mathbf{X}\_{t\_k} = \mathbf{x}\right] = \mathbb{E}\left[V(\mathbf{X}\_{t\_1}, t\_{k+1}) \mid \mathbf{X}\_{t\_0} = \mathbf{x}\right].$$

Thus, letting *Pt*<sup>1</sup> *v*(*x*) = **E** � *v*(*Xt*<sup>1</sup> ) | *Xt*<sup>0</sup> = *x* � be the expectation operator, from the backward induction (44) we get

$$\begin{cases} v\_M(\mathbf{x}) = \mathbf{1}\_{\{\mathbf{x} > l\}} \cdot K(1 - e^{\mathbf{x}})^+, \\ v\_k(\mathbf{x}) = \mathbf{1}\_{\{\mathbf{x} > l\}} P\_{l\_1} v\_{k+1}(\mathbf{x}), \quad j = M - 1, \dots, 1, \\\ v\_0(\mathbf{x}) = P\_{l\_1} v\_1(\mathbf{x}), \end{cases} \tag{45}$$

where *vk*(*x*) = *er*Δ*<sup>t</sup> V*(*x*, *tk*). And hence, the problem now becomes to find the function *v*0(*x*) from the backward induction (45).

Note that 1{*x*>0} <sup>=</sup> <sup>1</sup> 2 � 1 + sgn(*x*) � for all *x* ∈ **R**. Using the Property P4 (Relation with Hilbert transform) we obtain

$$\mathcal{F}(\mathbf{1}\_{\{x>0\}} \cdot v)(u) = \frac{1}{2} (\mathcal{F}v(u) + \mathrm{i}\mathcal{H}(\mathcal{F}v)(u)).\tag{46}$$

Let *a* ∈ **R** and T*<sup>a</sup>* be the transform operator: T*av*(*x*) = *v*(*x* − *a*). Then, we have

$$1\_{\{\mathbf{x} > \mathbf{a}\}} = \mathcal{T}\_{\mathbf{a}} 1\_{\{\mathbf{x} > \mathbf{0}\}} = \frac{1}{2} \left( 1 + \mathcal{T}\_{\mathbf{a}} \mathbf{sgn}(\mathbf{x}) \right) \lambda$$

and hence,

22 Will-be-set-by-IN-TECH

the authority to decide to execute the option now or reserve it to the next time point. So we

can split the integral that defines *Vj*(*tk*) into three parts: [*a*, *x*<sup>∗</sup>

*<sup>k</sup>* ; *tk*) + *Cj*(*x*<sup>∗</sup>

*<sup>k</sup>* ; *tML*) + *Cj*(*x*<sup>∗</sup>

**4.3 The fast Hilbert transform approach for pricing barrier options**

*<sup>G</sup>* <sup>=</sup> <sup>1</sup>{*St*

formula, the price of this option is given by

*V*(*x*, *tk*) = *e*−*r*Δ*<sup>t</sup>*

*V*(*x*, *t*0) = *e*−*rT***E**

**E** � <sup>−</sup>*rT***E** � <sup>1</sup>{*Xt*

*<sup>V</sup>*(*x*, *tM*) = <sup>1</sup>{*x*>*l*}*K*(<sup>1</sup> <sup>−</sup> *<sup>e</sup>x*)+,

*V*(*x*, *t*0) = *e*

⎧ ⎪⎪⎨

⎪⎪⎩

1, . . . , *M*−1, we have

an efficient computational algorithm based on the fast Hilbert transform.

which can be computed recursively by the following backward induction:

<sup>1</sup>{*x*>*l*}**<sup>E</sup>** �

> � *e*−*r*Δ*<sup>t</sup>*

*V*(*Xtk*<sup>+</sup><sup>1</sup> , *tk*<sup>+</sup>1) | *Xtk* = *x*

Let **T** = {*tk* : *k* = 1, . . . , *M*} be the set of pre-specified monitored dates, where

*<sup>k</sup>* < *h*, which means the early-exercise point doesn't hit the up barrier. Thus, We have

*<sup>k</sup>* ≥ *h* , which means the early-exercise point hits the up barrier. Thus, the option

� *Cj*(*a*, *<sup>h</sup>*; *tk*) + *<sup>e</sup>*−*r*(*T*−*tk*−<sup>1</sup>)*Dj*(*h*, *<sup>b</sup>*), for a call option, *Cj*(*a*, *h*; *tML*) + *e*−*r*(*T*−*tk*−<sup>1</sup>)*Dj*(*h*, *b*), for a put option,

Ding *et al* (2011a) gave a detail algorithm, which summarizes the above procedure, for pricing an up-and-out Bermudan barrier option, and the corresponding numerical experiments.

Feng & Linetsky (2008) presented a new numerical method to price discretely monitored barrier options under exponential Lévy models. Their method involves the relation with Hilbert transform (Property P4) and the Sine expansion in Hardy spaces. They also gave

<sup>0</sup> = *<sup>t</sup>*<sup>0</sup> < *<sup>t</sup>*<sup>1</sup> < ··· < *tM* = *<sup>T</sup>* with <sup>Δ</sup>*<sup>t</sup>* = *tk* − *tk*−<sup>1</sup> = *<sup>T</sup>*/*M*.

where *K* is the strike price and 0 < *L* < *K* is the lower barrier. We also assume that the underlying asset price process is given by *St* = *KeXt* with *S*<sup>0</sup> = *S*, where *Xt* is a Lévy process started at *x* = log(*S*/*K*). Denoting *l* = log(*L*/*K*), then, with help of the risk-neutral valuation

<sup>1</sup>>*l*}1{*Xt*2>*l*} ··· <sup>1</sup>{*XtM* <sup>&</sup>gt;*l*}*K*(<sup>1</sup> <sup>−</sup> *<sup>e</sup>*

*V*(*Xtk*<sup>+</sup><sup>1</sup> , *tk*<sup>+</sup>1) | *Xtk* = *x*

*V*(*Xt*<sup>1</sup> , *t*1) | *Xt*<sup>0</sup> = *x*

Since each Lévy process is stationary and has independent increments, for each *k* =

� = **E** � �

*V*(*Xt*<sup>1</sup> , *tk*<sup>+</sup>1) | *Xt*<sup>0</sup> = *x*

� .

<sup>1</sup>>*L*}1{*St*2>*L*} ··· <sup>1</sup>{*StM* <sup>&</sup>gt;*L*}(*<sup>K</sup>* <sup>−</sup> *ST*)+,

We consider a European-type barrier put option whose payoff at maturity *T* is given by

*<sup>k</sup>* ], (*x*<sup>∗</sup>

*<sup>k</sup>* , *<sup>h</sup>*; *tML*) + *<sup>e</sup>*−*r*(*T*−*tk*−<sup>1</sup>)*Dj*(*h*, *<sup>b</sup>*), for a call option,

*<sup>k</sup>* , *<sup>h</sup>*; *tk*) + *<sup>e</sup>*−*r*(*T*−*tk*−<sup>1</sup>)*Dj*(*h*, *<sup>b</sup>*), for a put option.

*<sup>k</sup>* , *h*) and [*h*, *b*]. We have

*XtM* )<sup>+</sup> <sup>|</sup> *<sup>X</sup>*<sup>0</sup> <sup>=</sup> *<sup>x</sup>*

, *k* = *M*−1, . . . , 1,

� . � ,

(44)

Case 1: *x*∗

Case 2: *x*∗

*Vj*(*tk*) =

� *Cj*(*a*, *<sup>x</sup>*<sup>∗</sup>

*Vj*(*tk*) =

*Cj*(*a*, *x*<sup>∗</sup>

integral can be split into two parts: [*a*, *h*) and [*h*, *b*]:

$$\mathbf{1}\_{\{\mathbf{x} > a\}} \cdot v(\mathbf{x}) = \frac{1}{2} \left( v(\mathbf{x}) + v(\mathbf{x}) \cdot \mathcal{T}\_a \text{sgn}(\mathbf{x}) \right) = \frac{1}{2} \left( v(\mathbf{x}) + \mathcal{T}\_a (\text{sgn} \cdot \mathcal{T}\_{-a} v)(\mathbf{x}) \right).$$

Taking the Fourier transform to both sides of this equation, we have

$$\mathcal{F}(1\_{\{x>a\}} \cdot v)(u) = \frac{1}{2}\mathcal{F}v(u) + \frac{1}{2}\mathcal{F}\left(\mathcal{T}\_a(\text{sgn} \cdot \mathcal{T}\_{-a}v)\right)(u).$$

and using equation (46) we obtain

$$\mathcal{F}(\mathbf{1}\_{\{\mathbf{x}>\mathbf{a}\}} \cdot \mathbf{v})(\boldsymbol{u}) = \frac{1}{2} \mathcal{F}\boldsymbol{v}(\boldsymbol{u}) + \frac{1}{2} \mathrm{i}\boldsymbol{e}^{\mathrm{i}\boldsymbol{u}\boldsymbol{u}} \mathcal{H}(\boldsymbol{e}^{-\mathrm{i}\boldsymbol{u}\boldsymbol{y}} \mathcal{F}\boldsymbol{v}(\boldsymbol{y}))(\boldsymbol{u}).\tag{47}$$

On other hand, noting that, for each *k* = 1, . . . , *M*−1, the condition density of *Xtk*<sup>+</sup><sup>1</sup> given *Xtk* = *x* possesses the property:

$$f(y|\mathbf{x}) = f(y-\mathbf{x}), \quad \mathbf{x}, y \in \mathbb{R}\_2$$

where *f*(*y*) is the density of *Xt*<sup>1</sup> under the initial condition *Xt*<sup>0</sup> = *x*. The infinite integrals *Pt*<sup>1</sup> *vk*<sup>+</sup>1(*x*) in (45) becomes to

$$P\_{l\_1}v\_{k+1}(\mathbf{x}) = \int\_{-\infty}^{+\infty} v\_{k+1}(y)f(y-\mathbf{x})d\mathbf{y} = \int\_{-\infty}^{+\infty} v\_{k+1}(\mathbf{x}+z)f(z)dz.$$

Then, this integral can be rewritten as a convolution of *vk*<sup>+</sup>1(*x*+*z*) and the function *f*(−*x*), i.e.

$$P\_{t\_1}v\_{k+1}(\mathbf{x}) = f(-\mathbf{x}) \* v\_{k+1}(\mathbf{x}).$$

Note that supp*vk*(*x*) ⊂ (*l*, 0) for each *k* = *M*, . . . , 1 from the backward induction (45). We can take the Fourier transform in the both sides of (45). Applying the Property P3 (Convolution) we have

$$\mathcal{F}\left[P\_{l\_1}v\_k(\mathbf{x})\right](\boldsymbol{u}) = \mathcal{F}\left[f(-\mathbf{x}) \* v\_{k+1}(\mathbf{x})\right](\boldsymbol{u}) = \mathcal{F}\left[f(-\mathbf{x})\right](\boldsymbol{u}) \cdot \mathcal{F}\left[v\_{k+1}(\mathbf{x})\right](\boldsymbol{u}).$$

Dempster, M. A. & Hong, S. S. G. (2002). Spread option valuation and the fast Fourier transform, *Mathematical Finance–Bachelier Congress 2000*, Spring, pp. 203–220. Ding, D & U, S. C. (2010). An accurate and stable FFT-based method for pricing options under

Fourier Transform Methods for Option Pricing 289

Ding, D & U, S. C. (2011). Efficient option pricing methods based on Fourier series expansions,

Ding, D., Huang, H. & Zhao, J. (2011a). An efficient algorithm for Bermudan barrier option

Ding, D., Weng, Z. & Zhao, J. (2011b). Pricing Bermudan barrier options via a fast and accurate FFT-based Method, to appear in the Proceedings of CiSE 2011, IEEE. Duffie, D., Pan, J. & Singleton, K. (2000). Transform analysis and asset pricing for affine

Dufresne, D., Garrido, J. & Morales, M. (2009). Fourier inversion formulas in option pricing

Eberlein, E., Glau, K., & Papapantoleon, A. (2010). Analysis of Fourier transform valuation formulas and application, *Applied Mathematical Finance*, Vol. 17, pp. 211–240. Fang, F. & Oosterlee, C. W. (2008). A novel pricing method for European option based on Fourier-cosine series expansions, *SIAM J. Sci. Comput.*, Vol. 31, pp. 826–848. Fang, F. & Oosterlee, C. W. (2009). Pricing early-exercise and discrete barrier options by

Fang, F. & Oosterlee, C. W. (2011). A Fourier-based valuation method for Bermudan and barrier options under Heston's model, *SIAM J. Financial Math.*, Vol. 2, pp. 439–463. Feng, L. & Linetsky, V. (2008). Pricing discretely monitored barrier options and defaultable

Gurland, J. (1948). Inversion formula for the distribution of ratios, *Ann. Math. Statist.*, Vol. 19,

Heston, S. (1993). A closed-form solution for options with stochastic volatility with applications to bond and currency options, *Rev. Financ. Stud.*, Vol. 6, pp. 327–343. Hurd, T. & Zhou, Z. (2010). A Fourier transform method for spread option pricing, *SIAM J.*

Lee, R. (2004). Option pricing by transform methods: Extensions, unification, and error

Lewis, A. (2001). A simple option formula for general jump-diffusion and other exponential

Lord, R., Fang, F., Bervoets, F. & Oosterlee, C. W. (2008). A fast and accurate FFT-based method

Kou, S. G. (2002). A jump-diffusion model for option pricing, *Management Science*, Vol. 48, pp.

Kwok, Y. K., Leung, K. S. & Wong, H. Y. (2010). Efficient options pricing using the fast Fourier

Madan, D., Carr, P. & Chang, E. (1998). The variance gamma process and option pricing,

Merton, R. (1973). Theory of rational option pricing, *Bell Journal of Economics and Management*

for pricing early-exercise options under Lévy processes. *SIAM J. Sci. Comput.*, Vol. 30,

control, *Journal of Computational Finance*, Vol. 7, pp. 51–86.

bonds in Lévy process models: a fast Hilbert transform approach. *Mathematical*

and insurance, *Methodol. Comput. Appl. Probab.*, Vol. 11, pp. 359–383.

Fourier-cosine series expansions, *Numer. Math.*, Vol. 114, pp. 27–62.

*Journal of Mathematical Research and Exposition*, Vol. 31, pp. 12–22.

pricing, Working paper, University of Macau.

jump-diffusions, *Econometrica*, Vol. 68, pp. 1343–1376.

Physics, pp. 741–746.

*Finance*, Vol. 18, pp. 337–384.

*Finan. Math.*, Vol. 1, pp. 142–157.

Lévy processes. http://www.optioncity.net

transform. http://ssrn.com/abstract=1534544

*European Financial Review*, Vol. 2, pp. 79–105.

*Science* , Vol. 4, pp. 141–183.

pp. 228–237.

pp. 1678–1705.

1086-1101.

exp-Levy processes, *The Proceedings of ISCM II - EPMSC XII*, American Institute of

Denote *ϕ*(*u*) is the characteristic function of *f*(*x*) on the complex plane **C**. Then, by a simple calculation, we have

$$\mathcal{F}\left[P\_{t\_1}v\_k(\mathbf{x})\right](\boldsymbol{u}) = \boldsymbol{\varrho}(-\boldsymbol{u}) \cdot \mathcal{F}v\_{k+1}(\boldsymbol{u}).\tag{48}$$

Thus, using the formulas (47) and (48) we obtain

⎧ ⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎩ *<sup>v</sup>*ˆ*M*(*u*) = *<sup>K</sup>*(1−*e*i*ul*) <sup>i</sup>*<sup>u</sup>* <sup>−</sup> *<sup>K</sup>*(1−*e*(1+i*u*)*<sup>l</sup>* ) <sup>1</sup>+i*<sup>u</sup>* , *v*ˆ*k*(*u*) = <sup>1</sup> <sup>2</sup> *<sup>ϕ</sup>*(−*u*)*v*ˆ*k*+1(*u*) + <sup>1</sup> <sup>2</sup> <sup>i</sup>*e*i*ul*H� *<sup>e</sup>*−i*yl <sup>ϕ</sup>*(−*y*)*v*ˆ*k*+1(*y*) � (*u*), *k* = *M*−1, . . . , 1, *v*ˆ0(*u*) = *ϕ*(−*u*)*v*ˆ1(*u*). (49)

Here *v*ˆ*k*(*u*) = F*vk*(*u*) for each *k*. Applying the truncated Sinc approximation, Feng & Linetsky (2008) obtained the discretization of *v*ˆ*k*(*u*):

$$\begin{split} \hat{\psi}\_{k}(nh) &= \frac{1}{2} \varphi(-nh) \hat{\psi}\_{k+1}(nh) \\ &+ \frac{\mathrm{i}e^{\mathrm{inhl}}}{2\pi} \sum\_{j=-N, j\neq n}^{N} e^{-\mathrm{i}jhl} \varphi(-jh) \hat{\psi}\_{k+1}(jh) \frac{1 - (-1)^{n-j}}{n-j}, \end{split} \tag{50}$$

for *n* = −*N*,..., *N* and *k* = *M*−1, . . . , 1, where *N* is a positive integer and *h* is the discretization step size. Then, the function *v*0(*x*) can be computed by the discretised inversion Fourier transform:

$$v\_0(\mathbf{x}) = \frac{1}{2\pi} \sum\_{j=-N}^{N} e^{-\mathbf{i}j\hbar\mathbf{x}} \varphi(-j\hbar) \vartheta\_1(j\hbar)\hbar.$$

Furthermore, Feng & Linetsky (2008) showed that the computation (50) involves a Toeplitz matrix-vector multiplication, which can be accomplished in *O*(*N* log2 *N*) operations using the FFT technique. They also referred the corresponding algorithm of computing the discrete Hilbert transform via the FFT as the *Fast Hilbert Transform*.

#### **5. References**


24 Will-be-set-by-IN-TECH

Denote *ϕ*(*u*) is the characteristic function of *f*(*x*) on the complex plane **C**. Then, by a simple

) <sup>1</sup>+i*<sup>u</sup>* ,

Here *v*ˆ*k*(*u*) = F*vk*(*u*) for each *k*. Applying the truncated Sinc approximation, Feng & Linetsky

for *n* = −*N*,..., *N* and *k* = *M*−1, . . . , 1, where *N* is a positive integer and *h* is the discretization step size. Then, the function *v*0(*x*) can be computed by the discretised inversion

Furthermore, Feng & Linetsky (2008) showed that the computation (50) involves a Toeplitz matrix-vector multiplication, which can be accomplished in *O*(*N* log2 *N*) operations using the FFT technique. They also referred the corresponding algorithm of computing the discrete

Barndorff-Nielsen, O. (1997). Normal inverse Gaussian distributions and stochastic volatility

Bates, D. S. (1996). Jumps and stochastic volatility: exchange rate processes implicit in

Black, F. & Scholes, M. (1973). The pricing of options and corporate liabilities, *Journal of Political*

Borak, S., Detlefsen, K. & Härdle, W. (2005). FFT based option pricing, SFB 649 Discussion

Carr, P., Geman, H., Madan, D. & Yor, M. (2002). The fine structure of asset returns: An

Carr, P. & Madan, D. (1999). Option valuation using the fast Fourier transform, *Journal of*

Carr, P. & Wu, L. (2003). The finite moment log stable process and option pricing, *Journal of*

Cont, R. & Tankov, P. (2003). *Financial Modelling with Jump Processes*, Chapman & Hall/ CRC

<sup>2</sup> <sup>i</sup>*e*i*ul*H�

<sup>−</sup>i*jhl <sup>ϕ</sup>*(−*jh*)*v*ˆ*k*+1(*jh*)

<sup>−</sup>i*jhxϕ*(−*jh*)*v*ˆ1(*jh*)*h*.

(*u*) = *ϕ*(−*u*) · F*vk*<sup>+</sup>1(*u*). (48)

� (*u*),

(49)

*k* = *M*−1, . . . , 1,

<sup>1</sup> <sup>−</sup> (−1)*n*−*<sup>j</sup>*

*<sup>n</sup>* <sup>−</sup> *<sup>j</sup>* , (50)

*<sup>e</sup>*−i*yl <sup>ϕ</sup>*(−*y*)*v*ˆ*k*+1(*y*)

F�

Thus, using the formulas (47) and (48) we obtain

*<sup>v</sup>*ˆ*M*(*u*) = *<sup>K</sup>*(1−*e*i*ul*)

*v*ˆ0(*u*) = *ϕ*(−*u*)*v*ˆ1(*u*).

+ i*e*i*nhl* 2*π*

*v*ˆ*k*(*u*) = <sup>1</sup>

(2008) obtained the discretization of *v*ˆ*k*(*u*):

*<sup>v</sup>*ˆ*k*(*nh*) = <sup>1</sup>

*Pt*<sup>1</sup> *vk*(*x*) �

<sup>i</sup>*<sup>u</sup>* <sup>−</sup> *<sup>K</sup>*(1−*e*(1+i*u*)*<sup>l</sup>*

<sup>2</sup> *<sup>ϕ</sup>*(−*u*)*v*ˆ*k*+1(*u*) + <sup>1</sup>

<sup>2</sup> *<sup>ϕ</sup>*(−*nh*)*v*ˆ*k*+1(*nh*)

*<sup>v</sup>*0(*x*) = <sup>1</sup>

Hilbert transform via the FFT as the *Fast Hilbert Transform*.

*Economy*, Vol. 81 (No. 3), pp. 637–654.

*Computational Finance*, Vol. 2 pp. 61–73.

*Finance*, Vol. 58, pp. 753–777.

Press.

paper 2005-011. http://sfb649.wiwi.hu-berlin.de

*N* ∑ *j*=−*N*,*j*�=*n*

2*π*

modelling, *Scandinavian Journal of Statistics*, Vol. 24, pp. 1–13.

Deutsche mark options, *Rev. Financ. Stud.*, Vol. 9, pp. 69–107.

empirical investigation. *Journal of Business*, Vol. 75, pp. 305–332.

*e*

*N* ∑ *j*=−*N e*

calculation, we have

Fourier transform:

**5. References**

⎧ ⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎩


**0**

**12**

Yi-Wen Liu

*Taiwan*

*National Tsing Hua University*

**Hilbert Transform and Applications**

Hilbert transform finds a companion function *y*(*t*) for a real function *x*(*t*) so that *z*(*t*) = *x*(*t*) + *iy*(*t*) can be analytically extended from the real line *t* ∈ R to upper half of the complex plane. In the field of signal processing, Hilbert transform can be computed in a few steps: First, calculate the Fourier transform of the given signal *x*(*t*). Second, reject the negative frequencies. Finally, calculate the inverse Fourier transform, and the result will be a complex-valued signal where the real and the imaginary parts form a Hilbert-transform pair.

When *x*(*t*) is narrow-banded, |*z*(*t*)| can be regarded as a slow-varying envelope of *x*(*t*) while the phase derivative *∂t*[tan−1(*y*/*x*)] is an instantaneous frequency. Thus, Hilbert transform can be interpreted as a way to represent a narrow-band signal in terms of amplitude and frequency modulation. The transform is therefore useful for diverse purposes such as latency analysis in neuro-physiological signals (Recio-Spinoso et al., 2011; van Drongelen, 2007), design of bizarre stimuli for psychoacoustic experiments (Smith et al., 2002), speech data compression for communication (Potamianos & Maragos, 1994), regularization of convergence problems in multi-channel acoustic echo cancellation (Liu & Smith, 2002), and

The rest of this review chapter is organized as follows: Sec. 2 reviews the mathematical definition of Hilbert transform and various ways to calculate it. Secs. 3 and 4 review applications of Hilbert transform in two major areas: Signal processing and system identification. The chapter concludes with remarks on the historical development of Hilbert

The desire to construct the Hilbert transform stemmed from this simple quest: Given a real-valued function *f* : R→R, can we find an imaginary part *ig* such that *fc* = *f* + *ig* can be analytically extended? For example, if *f*(*x*) = cos(*x*), then by inspection we can find *g*(*x*) = sin(*x*) such that *fc*(*x*) = *f* + *ig* = exp(*ix*). This function can obviously be extended analytically to the entire complex plane by replacing the real variable *x* with the complex

Re{ *f*ext(*z*)}|*z*=*<sup>x</sup>* = *f*(*x*), (1)

signal processing for auditory prostheses (Nie et al., 2006).

**2. Mathematical foundations of Hilbert transform**

variable *z* in the expression; the result is *f*ext(*z*) = exp(*iz*) and we have

**1. Introduction**

transform in Sec. 6.

