Fourier Transform Methods for Option Pricing *

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.

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.  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.

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: The Fourier transform of g(x) is defined by where i = √ −1 is the imaginary unit. Given the function F g(u), the function g(x) can be recovered by the Fourier inversion formula: Sometime it may be more convenient to consider the generalized Fourier transform on the complex plane. Let a, b ∈ R with a < b. Assume Then, the generalized Fourier transform: 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: with a < w < b. Here, Im{·} denotes taking the imaginary part of argument.
The following properties of Fourier transform are useful in this chapter.
P2.Modulation: F e λx g(x) (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), which is defined by 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: , where ·, · is the inner product of two square-integrable function defined by

Fast Fourier transforms
Let N = 2 L be a power of 2, and let x =(x 0 , x 1 ,...,x N−1 ) T be a given N-dimensional vector, where v T presents the transpose of a vector v. The discrete Fourier transform of x is another N-dimensional vector Dx =(y 0 , y 1 ,...,y N−1 ) T is defined by Denote A N an N × N matrix whose (p, j) th entry is given by Then, the discrete Fourier transform of x is given by It is clear that the computation to find Dx requires N 2 steps. However, the fast Fourier technique (FFT) would require only 1 2 NL = 1 2 N log 2 N steps. The idea behind the FFT algorithm is to take advantage of the periodicity property of N roots of unity. For the given It is easy to verify that the vector Dx, which is defined in (4), now is given by Instead of performing the matrix-vector multiplication A N x, we now only need to perform two matrix-vector multiplications A M x ′ and A M x ′′ so that the number of operations is reduced from N 2 to 2(N/2) 2 = 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 O(N 2 ) to O(N log 2 N). The similar FFT algorithm also can be used to calculate the discrete Fourier inversion transform: N jp y j , p = 0, 1, . . . , N−1.

Characteristic functions
Let X be a random variable having the density function f (x). Then, the characteristic function of X is defined by P8. X and Y have the same distribution function if and only if they have the same characteristic function.

Exponential Lévy models
An adapted stochastic process X t with X 0 = 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 X t be a Lévy process, and let B 0 be the set family of all Borel sets U ⊂ R whose closureŪ does not contain 0. Set is a σ-finite measure on B 0 , which is called the Lévy measure of X t . Then, the characteristic function of Lévy process X t has the Lévy-Khintchine representation: where the characteristic exponent function: Here the triple (α, σ 2 , ν(dx)) is called the Lévy triple of X t . 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} e x ν(dx) < ∞, the exponential e X t is a martingale if and only if 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 X t : There are many models in the financing modeling literature simply correspond to different choices for σ 2 and ν(dx) to the Lévy process X t .
The cases that σ 2 = 0 and ν(dx) is a finite measure, i.e. X t is a Lévy jump diffusion: where B t is a standard Brownian motion, {ξ k } is a sequence of independent and identically distributed random variables with common density f ξ (x), and N t is a Poisson process of intensity λ, such that B t , {ξ k } and N t are mutually independent. In these models, the Lévy triple is given by 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 N t , 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 μ J and standard deviation σ J , so that the characteristic function of X t is given by • The second one belongs to Kou (2002). In this model, ξ k s are non-symmetric double exponentially distributed, and the characteristic function of X t is given by where η and κ are parameters.
• The VG model, which is proposed by Madan et al (1998): X t is a variance gamma process with parameters σ, κ and θ, and its characteristic function is given by: • The NIG model, which is presented by Barndorff-Nielsen (1997): X t 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): X t is a CGMY process whose characteristic function is given by where the parameters C, G and M are nonnegative, Y < 2, and Γ(x) is the gamma function.
• The Finite Moment Log Stable model, which is considered by Carr & Wu (2003): X t is a finite moment log stable process whose characteristic function is given by where ω, α and σ are parameters.

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 S t is defined by the system of stochastic differential equations (SDEs): with the initial values S 0 = s 0 and V 0 = v 0 , where B S t and B V t are two standard Brownian motions with correlation ρ (i.e. B S , B V t = ρt), κ is the mean reversion, θ is the long-run variance level, and σ is the volatility-of-volatility parameter. Although S t is no longer a Lévy process, under the risk neutral probability the characteristic function of the log-price process X t = log S t is known in closed form: where x 0 = log s 0 , r is the interest rate, and γ = σ 2 (u 2 + iu)+(κ − iρσu) 2 . • The Baste's stochastic volatility model [Bates (1996)]: In this model, the price process S t and the variance process V t are given by the system of SDEs: where B S t and B V t are two standard Brownian motions with correlation ρ, and Z t is a compound Poisson process with intensity λ and log-normal distribution of jump sizes such that if k is its jump size then Similarly to the Heston's model, the characteristic function of X t = log S t is known in closed form: where the characteristic function of the diffusion part: and the characteristic function of the jump part:

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) : 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: 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 0 and strike price K is given by where r is the interest rate, T is the maturity, k = log K, and f T (x) is the density function of X T = log S T . It is clear that the second integral is the probability: Π 2 = P(X t > k). By the Fourier inversion formula (8) we have where ϕ T (u) is the characteristic function of X T . By a change of probability measure, we can verify that the second integral is give by Thus, the price of the European call now is given by a Black-Scholes type formula: 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(S T ) is given by where p T (x) is the density function of S T . By the Parseval relation we have

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 S t = e rt+X t , under a risk-neutral probability. Let f T (x) be the density of X T . Consider the price of a call option: where k = ln K is the log strike price. Note that V T (k) → S 0 = 1ask →−∞, and the function V T (k) is not square-integrable. Thus, we cannot express the Fourier transform in strike in terms of the characteristic function ϕ T (u) of X T and then find a range of strikes by Fourier inversion.
To obtain a square-integrable function we consider the modified call price: for some α > 0, which is chosen to improve the integrability. Carr & Madan (1999) showed that a sufficient condition for square-integrability ofṼ(k) is given by Consider the Fourier transform ofṼ(k : α): By (10) and the definition of characteristic functions we have where ϕ T (u) is the characteristic function of X T . Using the Fourier inverse transform we get Note that ψ T (v : α) is odd in its imaginary part and even in its real part. Since V(k) is real we have Now this integral can be evaluated by the numerical approximation using the trapezoidal rule: to apply the FFT algorithm, we set N as a power of 2, and define the grid points: with the step sizes Δk and Δv, which satisfy the Nyquist relation: ΔkΔv = 2π/N. Then, the numerical approximation can expressed as 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(S T ) at maturity T, where S t = S 0 e rt+X t is the stock price process under the risk neutral probability. To proceed, we assume that X t is a Lévy process having the analytic characteristic function ϕ T (z), which is regular in the strip whereS X is the complex conjugate set of S X .
Let s = log S 0 denote the logarithm of current stock value. Then, the option price is given by Under the assumption we can compute the generalized Fourier transform of V(s) by the Parseval relation: Option prices can now be given by the generalized Fourier inversion formula (3): with a < w < b, for a range of initial values s.
For a European call option, the function h(x)=( e rT+x − e k ) + , is Fourier integrable in the region Im{z} > 1, where its generalized Fourier transform can be computed explicitly: Hence, the generalized Fourier transform of call option price takes the form: and hence, the option price is given by 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 h for various claims, for instance, the put option, the covered call or the cash-secured put, etc.

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 where K 1 and K 2 are strike prices [see, e.g. Kwok et al (2010)]. Define X it = log S it and k i = log K i , i = 1, 2, and let f T (x 1 , x 2 ) be the joint density of X 1T and X 2T under the risk neutral probability. Then, the characteristic function of X 1T and X 2T is defined by the following two-dimensional Fourier transform: Consider the price of the correlation option: Following the Carr and Madam method we consider the modified call price: for some parameters α 1 > 0 and α 2 > 0, which are chosen such that this modified price is square integrable for negative value of k 1 and k 2 . Then, by a direct integral the Fourier transform ofṼ(k 1 , k 2 ) is given by , v 2 − i(α 2 + 1) (α 1 + iv 1 )(α 1 + 1 + iv 1 )(α 2 + iv 2 )(α 2 + 1 + iv 2 ) .
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 1T , S 2T ) by a combination of rectangular strips. (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.

Hurd & Zhou
h(x 1 , And denote Γ(z) the complex gamma function defined by Then, the method of Hurd and Zhou mainly relies on the following Fourier representation of the payoff function h(x 1 , x 2 ): where ǫ 1 and ǫ 2 are any two given real numbers with ǫ 2 > 0 and ǫ 1 + ǫ 2 < −1, and The proof of this formula can be found in the appendix of Hurd & Zhou (2010). Let X 1t = log S 1t and X 2t = log S 2t with X 10 = x 1 and X 20 = x 2 , and let ϕ T (u 1 , u 2 ) be the characteristic function of X 1T − x 1 and X 2T − x 2 . Then, we have E e i(u 1 X 1T +u 2 X 2T ) = e i(u 1 x 1 +u 2 x 2 ) ϕ T (u 1 , u 2 ). Now, using the formula (16), the price of the spread option (14) is expressed as an explicit two-dimensional Fourier inversion transform: 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.

Fourier expansion methods
Fang & Oosterlee (2008)  We recall that the price S t of underlying asset follows an exponential Lévy model. Let X t = log(S t /K), where K is the strike price. Consider the price of European call in the form: where f T (x) is the probability density of X T under the risk-neutral probability. Let ϕ T (u) be the characteristic function of f T (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.,φ and to consider the Fourier-cosine expansion of where

Then, we get an approximation of f T (x) in (18) by
Now, substituting (19) into (16), we obtain an approximation of the option price: where Φ n (0, b) and Ψ n (0, b) are two integrals given by: for any [c, d] ⊂ [a, b], which are analytically given by 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 where c 1 , c 2 , and c 4 are the first, second, and fourth cumulates of X t = log(S t /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.

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 be the set of pre-specified monitored dates and the set of pre-specified exercise dates, respectively, before maturity T, where Consider a discrete American barrier option, which is monitored at every t k ∈ T and can be exercised at each t k ∈ T e , namely Bermudan barrier option, whose payoff is given by Here S t k is the price of the underlying asset at time t k ∈ T, H > K is the constant barrier and R 0 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 S t k hits the barrier level H at one time t k ∈ T, and it can also be exercised at any time t k ∈ T e before maturity T.
Denote V(S, t k ) the value of this Bermudan barrier option at time t k and the spot price S t k = S. With help of the risk-neutral valuation formula, this price process can be computed recursively by the following backward induction: in specialty, the initial price is given 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 where X t is a Lévy process and S 0 is the initial price. Let f (·|x) be the condition density of X t k+1 given X t k = x for t k ∈ T. Set g(x)= (S 0 e x − K) + , for a call option, (K − S 0 e x ) + , for a put option.
Then the backward induction (23) can be rewritten by and v where v(x, t k )=V(S 0 e x , t k ) for any t k ∈ T, and h = log(H/S 0 ).
Since each Lévy process is stationary and has independent increments, the condition density f (·|x) possesses the property: where f (y) is the density of X t 1 under the initial condition X t 0 = x. Applying this property to infinite integrals c(x, t k ) in (24) it becomes to for any t k ∈ T. Then, this integral can be rewritten as a convolution of v(x+z, t k+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 (Convolution), the integral becomes 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Δt F e αx c(x, t k ) (u)=ϕ −(u−iα) ·F e αx v(x, t k+1 ) (u).
Thus, taking the inverse Fourier transform we obtain Denote where ϕ T (u) is the characteristic function of X T , 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 in (26): where Δu = 2π/κ. It is clear that these grids satisfy the Nyquist relation: ΔuΔy = 2π/N. Now, for each t k ∈ 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 ω n e iu j y n +αy n v(y n , t k+1 ) where the weights ω n are chosen as ω 0 = ω N−1 = 1 2 and ω n = 1 for n = 1, . . . , N−2. Noting that u 0 = − 1 2 NΔu and Δx = Δy = 2π/(NΔu), we have Now, we can employ the FFT algorithm to calculate the summations in the right side of (27). Once the integral c(x p , t k ) is computed, we can determine the early-exercise price S * t k , t k ∈ T e , by the procedure: for each t k ∈ T e , locate j k such that or, In the case (28) set x * (t k )=x j k , and in the case (29) set x * (t k )= 1 2 (x j k + x j k+1 ). Then the early-exercise price at every t k ∈ T e is given by S * t k = S 0 e x * (t k ) . 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.

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 S t follows an exponential Lévy model, and T and T e 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 where H > K is the constant barrier and R 0 is the contractual rebate. Denote V(S, t k ), t k ∈ T e , the value of this Bermudan barrier option at time t k and the spot price S t k = 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 Here S is the spot price of underlying asset, r > 0 is the interest rate. Let X t = log(S t /K) be the logarithm of the underlying asset price S t over the strike price K, and denote x = log(S/K) and h = log(H/K). Let f (·|x) be the condition density of X t k+1 given X t k = x for t k ∈ T. Set g(x)= K(e x − 1) + , for a call option, K(1 − e x ) + , for a put option.
Then the backward induction (23) and the price formula (30) can be rewritten by and v where v(x, t k )=V(Ke x , t k ) for any t k ∈ T.
We consider the infinite integrals c(x, t k ) 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 Fourier-cosine expansion on [a, b]: where w 0 = 1 2 and w j = 1 for all j = 1, 2, 3, . . .. We substitute this expansion into the integral (33) and then we can rewrite it bȳ where for each j = 1, 2, 3, . . ., and Since the Fourier-cosine expansion has a high accuracy with a few terms, we can truncate the infinity series (34) and approximatec(x, t k ) by leaving the first N terms, i.e., for any t k ∈ T, where w 0 = w N−1 = 1 2 and w j = 1 for j = 1,...,N−2. On the other hand, we can represent each integral (36) by approximation as the following Let ϕ(u; x) be the characteristic function of f (·|x). Then, we can approximate each F j (x) by And so, we get the further numerical approximations of the integralsĉ(x, t k ) in (37) bỹ In special, the approximation of initial price v(x, t 0 ) in (32) Meanwhile, since the characteristic function φ(u; x) possesses the property: where ϕ(u)=ϕ(u;0), the approximations (38) and (39) can be simplified tõ and the initial price of optioñ v(x, t 0 )=e −rΔt In order to use this approximate formulation we still need to compute the integrals V j (t k ). For convenience we introduce the following notions: for any a ≤ x 1 ≤ x 2 ≤ b, wherec(x, t ML )=g(x), andc(x, t k ), t k ∈ T, are given in (40).
And next, we consider to calculate the integrals V j (t k ) for t k ∈ T \ T e . We have V j (t k )=C j (a, h; t k )+e −r(T−t k−1 ) D j (h, b), j = 0, 1, . . . , N−1.
Since the integral D j (x 1 , x 2 ) has the analytic representation (45), we only need to calculate the integral C j (x 1 , x 2 ; t k ). Fang & Oosterlee (2009) gave an efficient numerical algorithm which approximates C j (x 1 , x 2 ; t k ) by using FFT method with the operation cost O(N log 2 (N)).
Finally, we consider to calculate the integrals V k (t k ) for t k ∈ T e . It is clear that we should find the valueṽ(x, t k ) in the last equation in (31), or equivalently, to determine the early-exercise point x * k at each time t k , which is the point where the continuation value is equal to the payoff, i.e.,c( Then, the problem becomes to find the root x * k of each equation h k (y)=0. Note that the functionc(y, t k ), 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 * k . Here if x * k 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 * k , t k ∈ T e , we have two different cases for an up-and-out barrier option: Case 1: x * k < 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 V j (t k ) into three parts: [a, x * k ], (x * k , h) and [h, b]. We have , for a put option.
Case 2: x * k ≥ 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]: , for a call option, , 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.

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 = {t k : k = 1, . . . , M} be the set of pre-specified monitored dates, where We consider a European-type barrier put option whose payoff at maturity T is given by 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 S t = Ke X t with S 0 = S, where X t 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 which can be computed recursively by the following backward induction: Since each Lévy process is stationary and has independent increments, for each k = 1, . . . , M−1, we have Thus, letting P t 1 v(x)=E v(X t 1 ) | X t 0 = x be the expectation operator, from the backward induction (44) we get where v k (x)=e rΔt V(x, t k ). And hence, the problem now becomes to find the function v 0 (x) from the backward induction (45).
Note that 1 {x>0} = 1 2 1 + sgn(x) for all x ∈ R. Using the Property P4 (Relation with Hilbert transform) we obtain Let a ∈ R and T a be the transform operator: T a v(x)=v(x − a). Then, we have 1 {x>a} = T a 1 {x>0} = 1 2 1 + T a sgn(x) , and hence, Taking the Fourier transform to both sides of this equation, we have F 1 {x>a} · v (u)= 1 2 F v(u)+ 1 2 F T a sgn ·T −a v (u). and using equation (46) we obtain On other hand, noting that, for each k = 1, . . . , M−1, the condition density of X t k+1 given X t k = x possesses the property: where f (y) is the density of X t 1 under the initial condition X t 0 = x. The infinite integrals P t 1 v k+1 (x) in (45) becomes to Then, this integral can be rewritten as a convolution of v k+1 (x+z) and the function f (−x), i.e.
Note that suppv k (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 Denote ϕ(u) is the characteristic function of f (x) on the complex plane C. Then, by a simple calculation, we have F P t 1 v k (x) (u)=ϕ(−u) ·Fv k+1 (u).
Herev k (u)=F v k (u) for each k. Applying the truncated Sinc approximation, Feng & Linetsky (2008) obtained the discretization ofv k (u): 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: Furthermore, Feng & Linetsky (2008) showed that the computation (50) involves a Toeplitz matrix-vector multiplication, which can be accomplished in O(N log 2 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. The book focuses on Fourier transform applications in electromagnetic field and microwave, medical applications, error control coding, methods for option pricing, and Helbert transform application. It is hoped that this book will provide the background, reference and incentive to encourage further research and results in these fields as well as provide tools for practical applications. It provides an applications-oriented analysis written primarily for electrical engineers, control engineers, signal processing engineers, medical researchers, and the academic researchers. In addition the graduate students will also find it useful as a reference for their research activities.