A radial basis collocation method for pricing American options under regime-switching jump-diffusion models

Applied Numerical Mathematics 65 (2013) 79–90

Contents lists available at SciVerse ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

A radial basis collocation method for pricing American options under regime-switching jump-diffusion models Ali Foroush Bastani ∗ , Zaniar Ahmadi, Davood Damircheli Department of Mathematics, Institute for Advanced Studies in Basic Sciences, P.O. Box 45195-1159, Zanjan, Iran

a r t i c l e

i n f o

Article history: Received 11 September 2011 Received in revised form 18 September 2012 Accepted 21 October 2012 Available online 3 December 2012 Keywords: Regime switching Lévy models Jump-diffusion process Radial basis functions American option pricing Coupled partial integro-differential equations (PIDEs)

a b s t r a c t The Markovian regime-switching paradigm has become one of the prevailing models in mathematical finance. It is now widely known that under the regime-switching model, the market is incomplete and so the option valuation problem in this framework will be a challenging task of considerable importance for market practitioners and academia. Our concern here is to solve the pricing problem for American options in a Markov-modulated jump-diffusion model, based on a meshfree approach using radial basis functions. In this respect, we solve a set of coupled partial integro-differential equations with the free boundary feature by expanding the solution vector in terms of radial basis functions and then collocating the resulting system of equations at some pre-specified points. This method exhibits a superlinear order of convergence in space and a linear order in time and also has an acceptable speed in comparison with some existing methods. We will compare our results with some recently proposed approaches. © 2012 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction Recent years have witnessed the widespread success of a series of asset price models, widely known under the heading of regime-switching models (see [5], [25] and [14] and the many references cited therein). These models are able to describe the behavior of a number of observed financial time series in a parsimonious, operational, tractable and precise way. A straightforward approach to introduce “regime change” into the model is to assume that the dynamics of the underlying financial asset is governed by a Lévy process whose parameters are dependent on the states of a continuous-time hidden Markov chain with states describing the regimes in which the financial or economic system operates. For a comprehensive overview of this subject we refer the reader to [6] and the references therein. Thanks to the realistic assumptions made in the formulation of these models and efficient calibration methods available for them, we are able to replicate some important market features like volatility smiles, skewness and term structure premiums by employing these models. It is now widely known that under the regime-switching model, the market is incomplete and so the option valuation problem in this framework will be a challenging task of considerable importance for market practitioners and academia. Among the many researchers that have addressed the option pricing problem under the regime-switching model, we must mention the following: Naik [26] provides an elegant analytic solution for pricing European-style options under a two-state regime-switching diffusion model. Elliott et al. [7] study the pricing of European options under Markov-modulated geometric Brownian motion using regime-switching random Esscher transform method. Elliott et al. [8] consider the pricing of European and American options under a Markov-modulated jump-diffusion model using generalized Esscher transform method and derive a set of coupled partial integro-differential equations (PIDEs) describing the option prices in different regimes.

*

Corresponding author. E-mail address: [email protected] (A. Foroush Bastani).

0168-9274/$36.00 © 2012 IMACS. Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.apnum.2012.10.005

80

A. Foroush Bastani et al. / Applied Numerical Mathematics 65 (2013) 79–90

Siu [30] considers the option pricing problem with a game theoretical approach and shows that the results are consistent with the ones obtained by the regime-switching Esscher transform method of Elliott and his co-workers. Boyarchenko and Levendorskiˇı [1] calculate the value of an American option using a generalized version of Carr’s randomization technique for regime-switching Lévy models. Yuen and Yang [32] present a trinomial tree method to price European options in a jumpdiffusion regime-switching model, extending the work of [26] to a multi-regime case. Jackson et al. [18,19] and Jaimungal and Surkov [20] devise a new Fourier space time-stepping (FST) algorithm which avoids the problems associated with finite difference and lattice-based methods by transforming the underlying system of PIDEs into Fourier space and then solving a linear system of ordinary differential equations (ODEs) in this space. Let us also mention the recent work of Huang et al. [17] who discuss three types of discretization schemes (explicit, implicit-direct control and implicit-penalty method) to approximate the resulting optimal control problem for pricing American options under a regime-switching model which they claim could be extended to the Markov-modulated jump-diffusion processes. The common drawback of some presently used methods which are based on finite differences or lattice approximations is that they are not in general easily implementable in higher dimensions and multi-regime cases and also suffer from low convergence rates and/or low accuracy of the results. On the other hand and for the promising FST approach of Jaimungal and Surkov [20] which works with the transformed problem in the Fourier space, alongside the gains in speed resulting from the application of FFT method, the difficulty in extending the method to multi-asset cases stimulates the researchers to search for other alternative approaches. In recent years, the need for reliable and efficient methods to solve partial differential equations (PDEs) arising from engineering and scientific applications, which are in general high dimensional and defined on complex geometries, has resulted in the development of new numerical ideas and powerful discretization schemes, generally known as the “meshfree” or “meshless” methods (see e.g. Fasshauer [9] for a good survey in this topic). As the name indicates, these methods fully eliminate the need for a base “mesh” or “grid” on which to approximate the solution and the only structure they require is the distribution of some “nodes” or “points” in the solution domain. Of particular importance in this category, is a collocationbased method which employs radial basis functions (RBFs) as a basis to expand the solution in terms of them [21]. The ease of implementation in more than two dimensions, high order of accuracy and flexibility of this method, have converted it to a favorite approach to solve many different and hard-to-solve problems specially in financial option pricing applications. As a representative list of works in this field, we mention the following: Hon and Mao [16] apply the global RBFs for spatial collocation and solve the European and American option pricing problems, when the underlying asset satisfies a Black– Scholes stochastic differential equation. Hon [15] develops a numerical method based on RBFs and quasi-interpolation to solve the Black–Scholes PDE for valuing American options. Fasshauer et al. [10] consider a meshfree approximation scheme for multi-asset American options in the context of penalty reformulation of the problem. Goto et al. [13] use RBF approximation to solve the Black–Scholes PDE in pricing European, barrier and Asian options. Pettersson et al. [27] propose an improved RBF-based method for pricing European basket options which is second-order accurate in time and spectrally accurate in space. Larsson et al. [24] use RBFs in conjunction with the generalized Fourier transform method for multidimensional European options. Chan and Hubbert [4] solve the PIDE arising from one-dimensional jump-diffusion processes (considering both Merton and Kou models) to price European vanilla calls/puts on dividend-paying stocks and also American puts. Golbabai et al. [12] use the front-fixing transformation as well as the RBF expansion to arrive at an implicit nonlinear system of first-order equations which is then solved by a predictor–corrector method. Recently, Saib et al. [29] have introduced a new RBF-based approach to solve American option pricing problem under the Merton jump-diffusion model which incorporates a differential quadrature to efficiently handle the boundary conditions. They solve the semi-discrete equations obtained after their approximation using an exponential time integration scheme and provide several numerical tests showing the superiority of their approach over the popular Crank–Nicolson method. In this paper, we present and analyze an RBF-based collocation scheme to solve a set of coupled PIDEs which describe the prices of an American put option in different states of the underlying hidden Markov chain. Working with the operator formulation of the problem, we first discretize the time variable by the backward Euler method and then use RBF collocation on the resulting space-only PDE to arrive at a set of algebraic equations, the solution of which gives us the option price vector in a continuous form. This makes our approach conceptually and methodologically different from the one presented in [17] which tries to solve a set of coupled partial differential equations (PDEs) variational inequalities (VIs) by a variety of implicit and explicit finite difference methods. Among other advantages of this method which will be discussed in Section 4, the distinctive feature of this approach is the ease of implementation in a meshfree framework to solve the pricing problem, when the number of dimensions and also the number of regimes is arbitrarily large. The remainder of this paper is organized as follows: In Section 2, we introduce the regime-switching Lévy models of interest. Section 3 contains the details of our radial basis function (RBF) methodology. In Section 4, we consider two option pricing problems as test cases to illustrate the efficiency of our proposed method, one already treated in the literature (see Jackson et al. [18,19]) and the other being a new higher dimensional one. Finally, and in Section 5, we conclude the paper by a discussion of ongoing and future lines of research. 2. Regime-switching Lévy processes Markov chains are frequently used for capturing random shifts between different states. In this section, we review the most important definitions from continuous-time Markov chains, Lévy processes with regime-switching (or

A. Foroush Bastani et al. / Applied Numerical Mathematics 65 (2013) 79–90

81

Markov-modulated) parameters and also option pricing in this framework (see Chourdakis [6] for a comprehensive treatment). Let αt be a continuous-time hidden Markov chain taking values among H different states. Each state represents a particular regime and is labeled by an integer between 1 and H , hence the state-space of αt is given by H = {1, . . . , H }. Let the matrix Q = (q i j ) H × H denote the generator of αt . From the classical Markov chain theory (see e.g. Brémaud [2]), we know that all entries of Qare non-negative, except the main diagonal which are all non-positive. Furthermore, we also have the relationship q ii = − j =i q ji . Discretizing the continuous-time Markov chain by partitioning the time into steps of size t, we obtain the infinitesimal transition probabilities:



P kl =

1 + qll t , k = l, qkl t , otherwise.

See [5] for more details. As noted in the previous section, introduction of a Markov chain in our formulation will result in an incomplete market [8]. It is also well known that the absence of arbitrage opportunity implies the existence of an equivalent martingale measure (EMM). However, this martingale measure is not unique if the market is incomplete. One of the main approaches to pick an EMM – which will be of particular importance in option pricing – is the time-honored Esscher transform, first introduced to the insurance pricing literature by Gerber and Shiu [11]. The details of finding the corresponding EMM for the incompleteness arising from the Markov-modulated jump-diffusion processes can be found in Elliott et al. [8] and Siu et al. [31]. Let W t be a standard Brownian motion defined on a risk-neutral probability space (Ω, F , P) and assume that this process is independent of the Markov chain αt . We consider the following regime-switching exponential Lévy model for describing the underlying asset price dynamics:

S t = S 0 e Xt . The log-price process X t will be constructed in the following manner: Consider a collection of independent Lévy processes {Y ti }iH=1 indexed by i. The increments of the log-price process will switch between these H different Lévy processes, depending on the state αt : α

d X t = dY t t . Each Lévy process Y ti is assumed to have a Lévy–Itô decomposition of the form



dY ti = μti dt + σti dW t +

zN i (dz, dt ),

i = 1, 2, . . . , H ,

(2.1)

R i

in which μ is the drift and σ i is the diffusion coefficient of the i-th Lévy process. In this equation, N i (·, t ) is a Poisson random measure defined on Borel subsets of R with ν i (.) as its associated Lévy measure, describing the discontinuities. We now consider the pricing of an American put option written on the underlying asset { S t }t 0 with strike price K and maturity date T . To obtain an equation with constant coefficients for the price of this option in each regime, we switch to log-prices and let x = log( S t / S 0 ). Then the transformed option price V i (x, t ) at time 0  t  T and regime αt = i will satisfy the following equation, due to the risk-neutral pricing principle (see for example Karatzas and Shreve [22]):

 +    X t = x, αt = i . V i (x, t ) = sup E e −r (τ −t ) K − e X τ t τ  T

(2.2)

In the above equation, τ is a stopping time satisfying t  τ  T and E is the expectation operator with respect to equivalent martingale measure P. This is the optimal stopping formulation of the American option pricing problem. We must note here that in writing these equations and all subsequent ones, the parameters μi and σ i in (2.1) are taken regime-independent and constant and also the fixed interest rate r is absorbed into the constant μ term to simplify the presentation of the material. One can now show that V i (x, t ) for i = 1, 2, . . . , H satisfy the following system of free boundary value problems [23]:

⎧ H  ⎪ ∂Vi ⎪ i ⎪ − L V + qi j V j = 0, x > x¯ i (t ), ⎪ i ⎪ ∂t ⎪ ⎪ j =1 ⎪ ⎪ ⎪ ⎪ V i (x, t ) = K − e x , x  x¯ i (t ), ⎪ ⎪ ⎪   ⎪ x + ⎪ ⎨ V i (x, T ) = K − e , lim V i (x, t ) = 0, x→∞ ⎪ ⎪ ⎪ ⎪ ⎪ lim V i (x, t ) = K − e x¯ i (t ) , ⎪ ⎪ x→¯xi (t ) ⎪ ⎪ ⎪ ⎪ ∂Vi ⎪ ⎪ lim (x, t ) = −1, ⎪ ⎪ ⎪ ∂x x →¯ x ( t ) i ⎩ x¯ i ( T ) = K ,

i = 1, 2, . . . , H , i = 1, 2, . . . , H , i = 1, 2, . . . , H , i = 1, 2, . . . , H , i = 1, 2, . . . , H , i = 1, 2, . . . , H , i = 1, 2, . . . , H ,

(2.3)

82

A. Foroush Bastani et al. / Applied Numerical Mathematics 65 (2013) 79–90

in which x¯ i (t ) for i = 1, . . . , H denote the optimal exercise boundaries and Li is the infinitesimal generator of the i-th Lévy process of the form

     1 1 Li g (x) = r + λi g (x) − r − σ 2 − λi ξ ∂x g (x) − σ 2 ∂xx g (x) − g (x + z, t )ν i (dz). 2

2

(2.4)

R

In this equation, λi stands for the jump intensity at state i and







e z − 1 dF ( z)

ξ= R

for the function F which is the distribution of jumps sizes. This is a set of coupled partial integro-differential equations with H free boundaries due to the regime-switching feature introduced in the underlying asset model. The analytical solution of the above system of PIDEs is not available at hand and so the need for efficient numerical approaches seems a necessity. In the sequel, we introduce our approach to solve this set of equations. 3. Numerical meshfree method Our aim in this section is to use a collocation scheme based on radial basis functions (RBFs) to find an approximate solution for the set of Eqs. (2.3). By using the change of variables τ = T − t and applying the backward Euler scheme in time, we can use the collocation method in each time step to find a continuous approximation in the whole interval. It is obvious that V i (x, τ ) for i = 1, . . . , H satisfy the following set of coupled PIDEs in operator form:

 ∂Vi + Li V i − q i j V j = 0, ∂τ H

i = 1, 2, . . . , H ,

(3.1)

j =1

which is valid in the space–time domain [−∞, ∞] × [0, T ]. In order to numerically approximate the solution, let us truncate the x-domain into the subdomain J = [xmin , xmax ] and discretize the time variable by setting τ n = nτ for n = 0, 1, . . . , N τ , in which τ = NT and then define τ





V in (x) ≡ V i x, τ n ,

n = 0, 1, . . . , N τ .

Now by applying the backward Euler scheme on (3.1), we arrive at the following equation

V in (x) − V in−1 (x)

τ

+ Li V in (x) −

H 

qi j V nj (x) = 0

j =1

for i = 1, 2, . . . , H and n = 0, 1, . . . , N τ which by rearrangement gives

V in (x) + τ Li V in (x) − τ

H 

qi j V nj (x) = V in−1 (x).

(3.2)

j =1

Let xmin = x1 < x2 < · · · < x N x = xmax be a set of distinct nodes in the interval J , which are not necessarily equispaced. We can now represent V in (x) as a linear combination of radial basis functions of the form

V in (x) =

Nx 





ρi(,nj) φ x − x j  ≡

j =1

Nx 

ρi(,nj) φ j (x).

(3.3)

j =1

By substituting (3.3) in (3.2) we obtain: Nx 

ρi(,nj) φ j (x) + τ

j =1

Nx  j =1

ρi(,nj) Li φ j (x) − τ

H 

qi j

Nx 

ρ (jn,s) φs (x) = V in−1 (x).

(3.4)

s =1

j =1

We now let V in, j be the result of collocating V in (x) at the point x j and let the vectors Vn , n = 0, 1, . . . , N τ , be defined by



n n Vn = V 1n,1 , . . . , V 1n, N x , . . . , V H ,1 , . . . , V H , N x

T

.

As a result of this collocation procedure, we arrive at the following linear system of equations

A. Foroush Bastani et al. / Applied Numerical Mathematics 65 (2013) 79–90

⎛

83

(n) ⎞

⎛ A + τ ( B 1 − q A ) ⎞ ρ −τ q12 A ... −τ q1H A 1 11 ⎜ (n) ⎟ 2 − τ q A A +  τ ( B − q A ) . . . − τ q A 21 22 2H ⎜ ⎟ ⎜ ρ2 ⎟ n −1 ⎜ ⎟ ⎟⎜ .. .. ⎝ ⎠ ⎜ .. ⎟ = V . ... . ⎝ . ⎠ −τ q H1 A −τ q H2 A . . . A + τ ( B H − q H H A ) ρ (n)

(3.5)

H

in which A i j = φ j (xi ),

B ki, j

(n)

= L φ j (xi ) and ρk k

(n) (n) (n) = [ρk,1 , ρk,2 , . . . , ρk, N x ] T .

In order to compute the entries of the matrices B 1 , B 2 , . . . , B H , we need to estimate the integral part in (2.4) and for this purpose we use a Gaussian quadrature formula in a bounded interval of the form [ zmin , zmax ] (which will be specified for each given Lévy measure in Section 4) to arrive at



zmax 

V in (x +

i

z)ν (dz) ≈ λ

i

R

V in (x +

≈ λi

Nx 

ρi(,nj)

j =1

zmin Nx 

z) f ( z) dz = λ

i

ρi(,nj)

j =1

m 

zmax  φ j (x + z) f ( z) dz zmin

ωk φ j (x + zk ) f (zk ),

(3.6)

k =1

in which the ωk ’s are the Gaussian quadrature coefficients. Using this last approximation, we fill out the entries of the coefficient matrix of (3.5) and having at hand the vector Vn−1 from the initial conditions or previous time step, we solve the above system by a standard direct linear solver (using the MATLAB lu command) in a step by step time-marching fashion. The final price of the American option at time level n for i-th regime will be of the form:

⎧ ⎨ max{( K − e x )+ , A ρi(n) }, for x = {xmin , xmax }, n V i (x) = K − e x , for x = xmin , ⎩ 0, for x = xmax , in which the last two cases are included to account for the necessary boundary conditions. In this manner, the price of the option today in all H states could be approximated. Since the current state of the world is unknown to the trader, he or she must decide, exogenously, on a probability pk that the world is in state k. Once these H probabilities are determined, the trader’s price for the option will be k=1 p k V k . 3.1. Treatment of the early exercise boundary Everyone knows that the early exercise feature of the American option contracts makes the pricing of them a complicated and cumbersome task. Specially in our pricing problem and being in the i-th state of the world, we must determine an optimal exercise boundary S¯ i (t ) from which to decide whether or not to exercise the option at time t by checking whether S (t ) <= S¯ i (t ) or S (t ) > S¯ i (t ) respectively. In terms of the log-price variable, this is equivalent to finding x¯ i (τ ) which based on the fifth equation in (2.3) amounts to solve the functional equation





V i x¯ i (τ ), τ = K − e x¯ i (τ )

(3.7)

for x¯ i (τ ) and then setting S¯ i ( T − τ ) = e x¯ i (τ ) . At time level n, we can approximate x¯ i (τn ) by finding the unique root of the equation F in (x) = V in (x) − K + x = 0 with an iterative scheme like the Newton method of the form (m+1)

xi ,n

(m)

= x(i ,mn) −

F in (xi ,n )

( F in ) (x(i ,mn) )

,

in which the function F in and its derivative ( F in ) at any x can be easily computed by the following expressions:

F in (x) =

Nx 

ρi(,nj) φ j (x) − K + x,

j =1





F in (x) =

Nx 

ρi(,nj) φ j (x) + 1.

j =1

This advantage in explicit determination of the early exercise boundary gives this approach a greater flexibility in comparison with other competing ones like finite difference and finite element methods.

84

A. Foroush Bastani et al. / Applied Numerical Mathematics 65 (2013) 79–90

Table 1 Commonly used radial basis functions. Type of basis function

φ(r )

Piecewise smooth RBFs Linear Cubic Generalized Duchon spline (GDS) Wendland

r r3 r 2k log r, k ∈ N r 2w , w > 0 and w ∈ /N (1 − r )k+ p (r ), p a polynomial, k ∈ N

Matérn

Γ (w) r

Infinitely smooth RBFs Multiquadric (MQ) Inverse Multiquadric (IMQ)

√

21− w

w

K w (r ), w > 0

r2 + c2

√

1

r 2 +c 2 2 2 e −c r

Gaussian

4. Numerical results In this section, we present two examples to better illustrate the use of the switching Lévy approach and the proposed pricing methodology in concrete situations. These examples are concerned with American put options in three and five world states respectively. In the first example, we assume that the stock price follows a Merton jump-diffusion process with an intensity parameter governed by a three-state hidden Markov chain. In the second one, we consider the Kou jumpdiffusion model with jump intensities having a discrete five-state Markov dynamics. Test problem (1). In this example, we assume a three-regime economy in which the dynamics of the underlying stock price in the i-th regime obeys a Merton jump-diffusion process with the Lévy measure

ν i ( z) = √

λi 2πσ j

 exp −

( z − μ j )2 2σ j2

 ,

where the intensity vector is given by

λ = [0.3, 0.5, 0.7]T , the generator matrix is defined by

 Q =

 −0.8 0.6 0.2 , 0.2 −1 0.8 0.1 0.3 −0.4

and the a priori state probabilities are given to be

p = [0.2, 0.3, 0.5] T . Other useful data are provided in the following table: S

K

T

σ

r

σj

μj

100

100

1

0.15

0.05

0.45

− 0.5

For this problem, we use a uniform distribution of points in the interval [xmin , xmax ] = [−6, 6] for the collocation process and truncate the integration domain in (3.6) according to



zmax = −2σ j2 log(σ j

√

 12

2 π /2 )

+ μ j,

zmin = − zmax , with  = 10e−12. We must note here that using these two bounds forces the total truncation error to be uniformly bounded by  and the derivation of them is described in full detail in [3] and [4]. Among the possible choices for the radial basis function φ (where some of the most commonly used ones are summarized in Table 1), we choose to work with cubic RBFs, due to their independence of any shape parameter, the optimal selection of which is a complicated task and under ongoing research study (see e.g. [28] and the references therein). We compare our results with the FST method proposed by Jackson et al. [18] the results of which are summarized in Table 2. As you can see from this table, increasing the number of nodes in space and time will result in better results in comparison with the figure 14.2502 which is reported with the FST method. Fig. 1 also displays graphically the option

A. Foroush Bastani et al. / Applied Numerical Mathematics 65 (2013) 79–90

85

Table 2 Value of American put option in the three world states. Nx

Nτ

Value

Change

128 256 512 1024 2048

64 128 256 512 1024

14.3128 14.2842 14.2664 14.2575 14.2539

2.86 × 10−3 1.78 × 10−3 0.88 × 10−3 0.36 × 10−3

log2 -ratio

Time (s)

0.7 1 1.30

0.57 3.90 20.62 218.64 1683.82

Fig. 1. American put option price curve at t = 0 for test problem (1).

price V i vs. the stock price S in each regime for the range [75, 200] computed at time t = 0. In Fig. 2, we have plotted option price surfaces that are obtained with this method for three different regimes. Based on the methodology described in Subsection 3.1, we have also computed the optimal exercise boundary in each regime and have plotted the results in Fig. 3. In order to illustrate the rate of decay of numerical errors when discretization parameters τ and x tend to zero, we have prepared Fig. 4 in which the vertical axis in both subfigures is the log relative error and is calculated by

 RBF  V − V FST  . V FST

Log Relative Error = log

It is evident from this figure that the rate of convergence of the RBF collocation scheme is linear in time and superlinear in space. To have a better sense of the convergence behavior of the proposed method with respect to the space variable, we have included a column entitled “log2 -ratio” in Table 2. The numbers in this column report the convergence power of the algorithm by measuring the price on successively finer grids. Assuming that

V approx = V exact + C (x) p , we can estimate the exponent p via the relation

p ≈ log2

V approx (x) − V approx (x/2) V approx (x/2) − V approx (x/4)

,

which converges to 1.3 as we refine the space grid. One may argue here that in spite of the maximal order of convergence for the cubic RBF interpolation which is 4, why we don’t observe the rate of convergence 2 (resulting from the second-order differentiation operator) in space? The reason could be attributed to the lack of regularity of the problem due to the free boundary feature which results in deteriorated convergence rate and other numerical difficulties in the solution process which needs further research to be resolved efficiently. Test problem (2). In this test problem, we assume that the stock price process follows the Kou jump-diffusion model where the jumps arrive at Poisson times and are distributed according to the law





ν i (z) = λi p η1 e−η1 z 1z0 + (1 − p )η2 eη2 z 1z<0 .

86

A. Foroush Bastani et al. / Applied Numerical Mathematics 65 (2013) 79–90

Fig. 2. The option price surface for test problem (1): (a) λ1 = 0.3, (b) λ2 = 0.5 and (c) λ3 = 0.7.

Fig. 3. Exercise boundary of American put option for test problem (1).

A. Foroush Bastani et al. / Applied Numerical Mathematics 65 (2013) 79–90

87

Fig. 4. Numerical accuracy in test problem (1): (a) Influence of the number of time steps N τ for fixed x = 0.01; (b) Influence of the number of space steps N x for fixed τ = 0.001. Table 3 Value of American put option for different intensity regimes and discretization parameters. RBF

FST

Nx

Nτ

λ1 = 0.1

λ2 = 0.3

λ3 = 0.5

λ4 = 0.7

λ5 = 0.9

256 512 1024

128 256 512

9.9053 9.8697 9.8565

10.4368 10.4020 10.3884

10.9731 10.9390 10.9252

11.5111 11.4772 11.4632

12.0442 12.0104 11.9962

8198

512

9.8505

10.3818

10.9182

11.4561

11.9890

We assume that our five-state Markov chain has a generator of the form:

⎛

−1 ⎜ 0.25 ⎜ Q = ⎜ 0.25 ⎝ 0.25 0.25

0.25 −1 0.25 0.25 0.25

0.25 0.25 −1 0.25 0.25

0.25 0.25 0.25 −1 0.25

⎞

0.25 0.25 ⎟ ⎟ 0.25 ⎟ ⎠ 0.25 −1

and that the economy switches between different jump intensities described by the vector

λ = [0.1, 0.3, 0.5, 0.7, 0.9]T . In this case, we suppose that the market could be in any of the five regimes with equal probability. Other corresponding information is given in the following table: S

K

T

σ

r

p

η1

η2

100

100

0.25

0 .5

0.05

0 .5

3

2

We use a uniform distribution of points in the interval [xmin , xmax ] = [−6, 6] as collocation points and use the following bounds for the truncation process in (3.6):

zmax = log( / p )/(1 − η1 ), zmin = − log





 /(1 − p ) /(1 − η2 ),

where we use the value of  = 10e−12. We refer the reader to [4] to see a full derivation of these bounds in order to obtain uniform truncation error bounds. Table 3 contains the option prices corresponding to each intensity regime reported for different values of N x and N τ . The convergence pattern of the figures in each column indicates that the method works well in this case. To show the performance of the method in the whole time interval, we have prepared Fig. 5 which corresponds to the option prices as a function of both S and t over the range [80, 200] × [0, 0.25]. Fig. 6 also illustrates the cross section of these surfaces with plane t = 0 for different regimes. The smooth and stable behavior of these curves and surfaces is another indicator of the usefulness of this approach in test cases with a higher number of regimes. Optimal exercise boundary of the American put option has also been shown in Fig. 7.

88

A. Foroush Bastani et al. / Applied Numerical Mathematics 65 (2013) 79–90

Fig. 5. The option price surface for test problem (2): (a) λ1 = 0.1, (b) λ2 = 0.3, (c) λ3 = 0.5, (d) λ4 = 0.7 and (e) λ5 = 0.9.

5. Concluding remarks In this paper, a radial basis function collocation approach is introduced to price American options in a regime-switching Lévy context. After a brief review of the process of deriving the set of coupled PIDEs describing the prices in different regimes, we present the details of our methodology which consists of first discretizing in time (by backward Euler scheme) and then collocating in space (by RBFs). The two test problems which are studied in this paper demonstrate that this

A. Foroush Bastani et al. / Applied Numerical Mathematics 65 (2013) 79–90

89

Fig. 6. American put option price curve at t = 0 for test problem (2).

Fig. 7. Exercise boundary of American put option for test problem (2).

approach has a good performance in terms of accuracy, stability and ease of implementation. This method can easily be extended to price other path-dependent exotics, like barrier options and American-style Asian options in the regime-switching Lévy framework. These works would open a new line of research in studying these products, as there is no documented pricing methodology for them, at least to the best of the author’s knowledge. This method could also be easily employed when there are regime shifts in other parameters of the model like the interest rate and volatility. Acknowledgements The authors would like to thank the Editor and anonymous reviewers for their valuable comments and suggestions which were helpful in improving the paper. The first author gratefully acknowledges support by the Institute for Advanced Studies in Basic Sciences (IASBS) Research Council under grant No. G2011IASBS111. References [1] [2] [3] [4]

S. Boyarchenko, S. Levendorskiˇı, American options in regime-switching models, SIAM Journal on Control and Optimization 48 (2009) 1353–1376. P. Brémaud, Markov Chains: Gibbs Fields, Monte Carlo Simulation and Queues, Springer-Verlag, 1999. M. Briani, R. Natalini, G. Russo, Implicit–explicit numerical schemes for jump-diffusion processes, Calcolo 44 (2007) 33–57. R.T.L. Chan, S. Hubbert, A numerical study of radial basis function based methods for option pricing under one dimension jump-diffusion model, Applied Numerical Mathematics, submitted for publication, http://arxiv.org/abs/1011.5650. [5] K. Chourdakis, Continuous-time regime switching models and applications in estimating processes with stochastic volatility and jumps, Econometrica, submitted for publication, http://ssrn.com/abstract=358244, http://dx.doi.org/10.2139/ssrn.358244. [6] K. Chourdakis, Switching Lévy models in continuous-time: finite distributions on option pricing, preprint, http://ssrn.com/abstract=838924.

90

A. Foroush Bastani et al. / Applied Numerical Mathematics 65 (2013) 79–90

[7] R.J. Elliott, L. Chan, T.K. Siu, Option pricing and Esscher transform under regime switching, Annals of Finance 1 (2005) 423–432. [8] R.J. Elliott, T.K. Siu, L. Chan, J.W. Lau, Pricing options under a generalized Markov modulated jump diffusion model, Stochastic Analysis and Applications 25 (2007) 821–843. [9] G. Fasshauer, Meshfree Approximation Methods with MATLAB, World Scientific, 2007. [10] G.E. Fasshauer, A.Q.M. Khaliq, D.A. Voss, Using meshfree approximation for multi-asset American option problems, Journal of the Chinese Institute of Engineers 27 (2004) 563–571. [11] H.U. Gerber, E.S.W. Shiu, Option pricing by Esscher transforms, Transactions of Society of Actuaries 46 (1994) 99–191. [12] A. Golbabai, D. Ahmadian, M. Milev, Radial basis functions with application to finance: American put option under jump diffusion, Mathematical and Computer Modelling 55 (2011) 1354–1362. [13] Y. Goto, Z. Fei, S. Kan, E. Kita, Options valuation by using radial basis function approximation, Engineering Analysis with Boundary Elements 31 (2007) 836–843. [14] J.D. Hamilton, Regime-switching models, in: New Palgrave Dictionary of Economics, 2nd ed., Palgrave Macmillan, 2008. [15] Y.C. Hon, A quasi-radial basis functions method for American option pricing, Computers and Mathematics with Applications 43 (2002) 513–524. [16] Y.C. Hon, X.Z. Mao, A radial basis function method for solving options pricing model, Financial Engineering 8 (1999) 31–49. [17] Y. Huang, P.A. Forsyth, G. Labahn, Methods for pricing American options under regime switching, SIAM Journal on Scientific Computing 33 (2011) 2144–2168. [18] K.R. Jackson, S. Jaimungal, V. Surkov, Option pricing with regime switching Lévy processes using Fourier space time stepping, in: Proceedings of the Fourth IASTED International Conference on Financial Engineering and Applications, 2007, pp. 92–97. [19] K.R. Jackson, S. Jaimungal, V. Surkov, Fourier space time-stepping for option pricing with Lévy models, The Journal of Computational Finance (2009) 1–29. [20] S. Jaimungal, V. Surkov, Stepping through Fourier space, Risk (2009) 78–83. [21] E.J. Kansa, Multiquadrics – a scattered data approximation scheme with applications to computational fluid dynamics, I: Surface approximations and partial derivative estimates, Computers and Mathematics with Applications 19 (1990) 127–145. [22] I. Karatzas, S.E. Shreve, Methods of Mathematical Finance, Springer-Verlag, 1998. [23] A.Q.M. Khaliq, R.H. Liu, New numerical scheme for pricing American option with regime-switching, International Journal of Theoretical and Applied Finance 12 (2009) 319–340. [24] E. Larsson, K. Ahlander, A. Hall, Multi-dimensional option pricing using radial basis functions and the generalized Fourier transform, Journal of Computational and Applied Mathematics 222 (2008) 175–192. [25] R.S. Mamon, R.J.E. Elliott, Hidden Markov Models in Finance, Springer, 2007. [26] V. Naik, Option valuation and hedging strategies with jumps in the volatility of assets returns, Journal of Finance 48 (1993) 1969–1984. [27] U. Pettersson, E. Larsson, G. Marcusson, J. Persson, Improved radial basis function methods for multi-dimensional option pricing, Journal of Computational and Applied Mathematics 222 (2008) 82–93. [28] C.M.C. Roque, A.J.M. Ferreira, Numerical experiments on optimal shape parameters for radial basis functions, Numerical Methods for Partial Differential Equations 26 (2010) 675–689. [29] A.E.F. Saib, D.Y. Tangman, M. Bhuruth, A new radial basis functions method for pricing American options under Merton’s jump-diffusion model, International Journal of Computer Mathematics (2012), http://dx.doi.org/10.1080/00207160.2012.690034. [30] T.K. Siu, A game theoretic approach to option valuation under Markovian regime-switching models, Insurance: Mathematics and Economics 42 (2008) 1146–1158. [31] T.K. Siu, J.W. Lau, H. Yang, Pricing participating products under a generalized jump-diffusion model, Journal of Applied Mathematics and Stochastic Analysis (2008), http://dx.doi.org/10.1155/2008/474623. [32] F.L. Yuen, H. Yang, Option pricing in a jump-diffusion model with regime switching, ASTIN Bulletin 39 (2009) 515–539.