On the use of boundary conditions for variational formulations arising in financial mathematics

Applied Mathematics and Computation 124 (2001) 197±214 www.elsevier.com/locate/amc

On the use of boundary conditions for variational formulations arising in ®nancial mathematics Michael D. Marcozzi a, Seungmook Choi b, C.S. Chen c,* a

Department of Mathematical Sciences, University of Nevada Las Vegas, Las Vegas, NV, USA b Department of Finance, University of Nevada Las Vegas, Las Vegas, NV, USA c Department of Mathematical Sciences, University of Nevada Las Vegas, Box 454020, 4505 Maryland Parkway, Las Vegas, NV 89154-4020, USA

Abstract The general intractability of derivative security valuation models to present techniques, both analytic and numerical, arguably remains one of the preeminant problem of mathematical ®nance. It is the focus of this paper to examine the applicability of a promising recent development, namely Radial Basis Functions (RBF), to the problem of option valuation. A Black±Scholes framework is considered for American and European options written on a one and two risky assets. The performance of RBF and Finite-Di€erencing algorithms are examined with respect to arti®cial boundary conditions, computational domain, domain decomposition, and mesh scaling. Ó 2001 Elsevier Science Inc. All rights reserved. Keywords: Black±Scholes equation; Radial basis function; Domain decomposition; Bermuda approximation; Arti®cial boundary conditions

1. Introduction Over the last two decades, increasingly sophisticated derivative securities and securities with option components have been introduced by business corporations and ®nancial markets. Unfortunately, closed form solutions for *

Corresponding author. E-mail address: [email protected] (C.S. Chen).

0096-3003/01/$ - see front matter Ó 2001 Elsevier Science Inc. All rights reserved. PII: S 0 0 9 6 - 3 0 0 3 ( 0 0 ) 0 0 0 8 7 - 4

198

M.D. Marcozzi et al. / Appl. Math. Comput. 124 (2001) 197±214

such complex valuations rarely exist, particularly when the holder may e€ect cash ¯ows. Accordingly, numerical approximation techniques have been employed. Lattice techniques, of course, remain ubiquitous in practice for single-asset valuations, primarily as a result of their pedagogical simplicity and ease of implementation. Generalization of the lattice framework to higher dimensions have been introduced in [6,19]. The so-called stochastic mesh techniques, which serve, in e€ect, to generalize the deterministic lattice framework have recently been introduced in [7]. Unlike lattice techniques, variational formulations have received relatively less attention in the ®nance literature. We remark that the speci®cation of boundary conditions associated with the pricing equation has to date limited the application of this approach; in particular for multi-shock applications. Indeed, only through the use of accuracy limiting, application speci®c, ad hoc boundary conditions have such computations been performed for a single-asset (cf. [27]) or recently for two-shocks (cf. [28]). The application of arti®cial boundary conditions have, to our knowledge, never been analysed to ascertain their impact on method convergence (e.g. it appears to be ill-conceived to re®ne a mesh when no further improvements in accuracy may be attained due to the approximation properties of the arti®cial boundary conditions). In this paper, we consider the prototype problem of option valuation in a Black±Scholes framework for both a single-asset and two-assets relative to a variational formulation. We remark that this is in no way a restriction on the analysis (cf. [5,16,21,23]). This paper makes the following contributions: (i) we solve the variational pricing equations without introducing any boundary conditions, (ii) using arti®cial boundary conditions as well as a ``no boundary condition'' format, we investigate the e€ect of the computational domain radius on asymptotic convergence, (iii) we present an easily implemented general temporal error control mechanism, and (iv) we introduce a promising new technology for solving variational problems, namely Radial Basis Functions (RBF) (cf. [18,24]). Additionally, we discuss (v) domain decomposition techniques which are essential to RBF and ®nite di€erencing mesh re®nements as well as (vi) ``mesh scaling'' which allows a mesh to be ``thinned'' asymptotically for optimal eciency. We remark that to our knowledge, the topics (i), (ii), (iv), and (vi) are completely new, while (iii) and (v) are new to the ®nance literature. The primary advantage of RBF over existing variational techniques such as ®nite-di€erencing and ®nite-elements is their potential relative eciency approximating reasonably high dimensional problems (cf. [13]). A second advantage lies in their apparent simplicity and ease of implementation. Finally, as RBF represent a ``global'' interpolation technique, our computations show that they appear to be signi®cantly more stable than the second-order ®nite di€erencing methods in the absence of boundary conditions when applied to the initial-value formulations of optimal stopping problems.

M.D. Marcozzi et al. / Appl. Math. Comput. 124 (2001) 197±214

199

The outline of this paper is as follows. In Section 2, we consider the approximation of the one-asset option problem by RBF in discrete linear complementary form. In Section 3 the e€ects of mesh size, time discretization, arti®cial boundary conditions, and computational domain are examined with respect to the Black±Scholes pricing formula for both RBF and ®nite di€erence approximations. Domain decomposition and mesh scaling algorithms are implemented for the European and American cases. The two-asset case is considered in Section 4 and a summary of these ®ndings is presented in Section 5. 2. The single-asset case We consider in this section the application of RBF to the approximation of vanilla options dependent upon a single risky asset (cf. [2,4,20]). The discretization is presented for both European and American cases for the deterministic variational formulations associated with the valuation problems. To this end, we suppose a Black±Scholes framework wherein for a ®xed risk-free rate of return r, the asset price St satis®es the stochastic di€erential equation dSt ˆ rSt dt ‡ rSt dBt ; where t represents time, r2 the volatility and dBt so-called Brownian motion with respect to the risk-neutral measure. 2.1. European put A European option is de®ned as a non-negative, adapted process fh…t†g0 6 t 6 T such that h…T † is the payo€ of the claim at expiry T. For the prototypical example of the European put (call) with exercise price K, one has h…T † ˆ …K ST †‡ (h…T † ˆ …ST K†‡ ). 1 The arbitrage-free price of the put at time t is then given by P …St ; t† ˆ e

r…T t†

E‰…K

‡

ST † j Ft Š;

where fFt g0 6 t 6 T is the ®ltration associated with the Brownian motion and E the expectation with respect to the risk-neutral measure. Formally, the put price P …S; t† may be represented as a solution of the Black±Scholes equation for S > 0, and t 2 …0; T †: oP 1 2 2 o2 P …S; t† rS …S; t† ot 2 oS 2 with the terminal data P …S; T † ˆ …K 1

rS

oP …S; t† ‡ rP …S; t† ˆ 0; oS

‡

S† :

We use the notation …a†‡ :ˆ maxfa; 0g throughout.

…2:1a†

…2:1b†

200

M.D. Marcozzi et al. / Appl. Math. Comput. 124 (2001) 197±214

In order to set (2.1a) into canonical form, we make the following non-dimensional substitutions (cf. [14,27]): Let 2 s S ˆ Kex ; t ˆ T ; P ˆ Ke …1=2†…j 1†x …1=4†…j‡1† s p…x; s†; 2 …1=2†r where jˆ

r : …1=2†r2

The non-dimensional put value then satis®es the di€usion equation op …x; s† os

o2 p …x; s† ˆ 0 ox2

…2:2a†

for all …x; s† 2 R  R‡ such that  ‡ p…x; 0† ˆ e…1=2†…j 1†x e…1=2†…j‡1†x

…2:2b†

for all x 2 R. For a given mesh fx1 ; x2 ; . . . ; xN g; we suppose p…x; s† ˆ

N X

aj …s†
/…rj …x††

…2:3a†

jˆ1

such that /…rj …x†† ˆ ‰rj2 ‡ e c 2 Š1=2 ;

rj …x† ˆ jx

xj j;

…2:3b†

where e c is a constant known as the shape parameter. For a uniform mesh with spacing h ˆ xi xi 1 we take e c ˆ 4 h; although e c may be optimized for particular problems (cf. [8,18]). Substituting (2.3a) and (2.3b) into (2.2a) and (2.2b), it follows that for all s 2 R‡ and i ˆ 1; 2; . . . ; N : " # 2 N N X X daj 1 …xi xj † …s† /…rij † aj …s† ˆ 0; …2:4a† /…rij † ds /3 …rij † jˆ1 jˆ1 such that N X

 aj …0†/…rij † ˆ e…1=2†…j

1†xi

e…1=2†…j‡1†xi

‡

;

…2:4b†

jˆ1

where rij :ˆ ri …xj †. We remark that (2.4a) constitutes a coupled system of ®rstorder ordinary di€erential equations in aj …s†. In particular, the matrix U ˆ ‰/…rij †Š is symmetric positive de®nite and so invertible. It follows then that the representation (2.3a), with aj …0† uniquely determined by (2.4b), satis®es (2.4a) for all s 2 R‡ and so serves as the RBF approximation of the solution of (2.2a),

M.D. Marcozzi et al. / Appl. Math. Comput. 124 (2001) 197±214

201

the non-dimensionalized Black±Scholes equation. This is a re¯ection of the collocation (2.3a) and the semigroup properties of solutions to (2.2a) (cf. [14]). To e€ect comparisons relative to uniform time discretizations, we let …m†

aj daj …s† ˆ ds

…m 1†

aj Ds

…2:5†

for m ˆ 1; 2; . . . ; M, where Ds > 0, …m†

aj

ˆ aj …m
Ds†

and M
Ds ˆ 12r2 T . Substituting (2.5) into (2.4a), we obtain for i ˆ 1; 2; . . . ; N : " # N N N X X 1 1 …xi xj †2 …m† X 1 …m† …m 1† /…rij †aj /…rij †aj 0ˆ a j 3 Ds /…r Ds / …rij † ij † jˆ1 jˆ1 jˆ1 or, in operator form   1 U L a…m† 0ˆ Ds

1 Ua…m Ds

1†

:ˆ A a…m†

a:

…2:6†

By the invertibility of U ˆ ‰/…rij †Š and the speci®cation of aj …0† through (2.4b), the system (2.6) may be solved successively for each 1 6 m 6 M. 2.2. American put An American option is de®ned as a non-negative, adapted process fh…t†g0 6 t 6 T , where h…t† is the payo€ of the claim if exercised at time t. For the prototypical example of the American put (call) with exercise price K, one has ‡ ‡ h…t† ˆ …K St † (h…t† ˆ …St K† ). The arbitrage-free price of the put at time t is then given by   P …St ; t† ˆ inf E e r…h t† …K Sh †‡ j Ft ; t6h6T

where again fFt g0 6 t 6 T is the ®ltration associated with the Brownian motion and E the expectation with respect to the risk-neutral measure. Formally, the put price P …S; t† may be represented as a solution of the variational inequality for S > 0, and t 2 …0; T †:  oP 1 2 2 o2 P oP …S; t† rS …S; t† min …S; t† rS 2 ot 2 oS oS …2:7a†  ‡ ‡ rP …S; t†; P …S; t† …K S† ˆ0 such that P …S; T † ˆ …K

S†‡ ;

S > 0:

…2:7b†

202

M.D. Marcozzi et al. / Appl. Math. Comput. 124 (2001) 197±214

Nondimensionalizing (2.7a) and (2.7b) as for the European case, we have then that  op o2 p …x; s† min …x; s†; p…x; s† os ox2 …2:8a†    …1=4†‰…j 1†2 ‡4jŠs …1=2†…j 1†x …1=2†…j‡1†x ‡ e e ˆ0 e for all …x; s† 2 R  R‡ such that  ‡ p…x; 0† ˆ e…1=2†…j 1†x e…1=2†…j‡1†x ;

x 2 R:

…2:8b†

The linear complementary form associated with (2.8a) and (2.8b) is then    op o2 p …x; s† …x; s† p…x; s† ej1 s ‰ej2 x ej3 x Š‡ ˆ 0 …2:9a† 2 os ox for all …x; s† 2 R  R‡ such that o2 p …x; s† P 0; ox2

op …x; s† os and

 p…x; 0† ˆ ej2 x

where j1 ˆ

1h …j 4

j x

e3

‡

;

i 1†2 ‡ 4j ;

p…x; s†

ej1 s ‰ej2 x

‡

ej3 x Š P 0

x 2 R; 1 j2 ˆ …j 2

…2:9b†

…2:9c†

1†;

1 j3 ˆ …j ‡ 1†: 2

For a given mesh fx1 ; x2 ; . . . ; xN g; we again suppose that p…x; s† ˆ

N X

aj …s†
/…rj …x††;

jˆ1

where

h i1=2 /…rj …x†† ˆ rj2 ‡ ec 2 ;

rj …x† ˆ jx

xj j;

in which case for all s 2 R‡ and i ˆ 1; 2; . . . ; N : ( " #) 2 N N X X daj 1 …xi xj † …s†/…rij † aj …s† /…rij † ds /3 …rij † jˆ1 jˆ1 ( ) N X j1 s j2 xi j3 xi ‡ aj …s†/…rij † e ‰e e Š ˆ0  jˆ1

…2:10a†

M.D. Marcozzi et al. / Appl. Math. Comput. 124 (2001) 197±214

such that N X daj …s†/…rij † ds jˆ1 N X

N X jˆ1

"

1 aj …s† /…rij †

ej1 s ‰ej2 xi

aj …s†/…rij †

2

…xi xj † /3 …rij †

203

# P 0;

…2:10b†

‡

ej3 xi Š P 0

…2:10c†

jˆ1

and N X

aj …0†/…rij † ˆ ‰ej2 x

‡

ej3 x Š ;

x 2 R:

…2:10d†

jˆ1

Discretizing (2.10a)±(2.10d) uniformly in time, we again let …m†

aj daj …s† ˆ ds

…m 1†

aj Ds

for Ds > 0 and m ˆ 1; 2; . . . ; M, where as before …m†

aj

ˆ aj …m
Ds†;

in which case we de®ne for i ˆ 1; 2; . . . ; N : " # 2 N N X X 1 1 …xi xj † …m† …m† /…rij †aj wi :ˆ aj 3 Ds /…r † / …r † ij ij jˆ1 jˆ1 or in operator form   1 U L a…m† w :ˆ Ds

1 Ua…m Ds

1†

:ˆ Aa…m†

a

N X 1 …m /…rij †aj Ds jˆ1

1†

…2:11a†

and zi :ˆ

N X jˆ1

…m†

/…rij †aj

ej1 …m
Ds† ‰ej2 xi

‡

ej3 xi Š

or z :ˆ Ua…m†

b:

…2:11b†

We may then de®ne the discrete linear complimentary problem …q; M† as follows; ®nd w; z satisfying wt
z ˆ 0;

w P 0;

z P 0;

…2:12a†

where w ˆ Mz ‡ q

…2:12b†

204

M.D. Marcozzi et al. / Appl. Math. Comput. 124 (2001) 197±214

such that M ˆ AU

1

and q ˆ AU 1 b

ˆ

1 I Ds

LU

1

a:

Due to the invertibility of U ˆ ‰/…rij †Š, it follows that there exists a unique solution to …q; M†, for all Ds suciently small (cf. [10,22,25]), and consequently the RBF discretization is uniquely solvable. We remark that a so-called Bermuda approximation to the American formulation would be to solve (2.10d) for a…0† , time step Dt solving (2.6) for a…1† , and then applying the constraint (2.10c). This procedure may then be continued inductively for m ˆ 2; 3; . . . ; M. It subsequently follows from Lemke's algorithm that the so-called Bermuda approximation represents a ``®rst'' complementary approximation to the American solution. We note also that upon letting w  0 we recover Eq. (2.6); the European approximation. 3. Computational results In this section, we examine the implementation of both the RBF approximation and a traditional second-order central-space, implicit ®rst-order time, ®nite di€erencing scheme (FD) for the single-asset case (cf. [16,27]). Computational asymptotics are provided for the European case. In particular, we examine the dependency of the numerical approximations on the choice of mesh size, computational domain, and arti®cial boundary conditions. We remark that such computations are presented in order to allow an assessment of the limitations and bene®ts of the two approaches. Formulations utilizing mesh scaling and domain decomposition are also considered. The e€ectiveness of active time step management is also considered. Numerical solutions for the American case are presented. All estimates are considered relative to the L2 -norm; that is, for any compact set (i.e. closed and bounded) G  R and real-valued function f 2 C…G†; Z 2 2 kf kL2 …G† :ˆ ‰f …x†Š dx: G

The integral is approximated by the trapezoidal rule. All computations are performed relative to the non-dimensionalized formulations and symmetric computational and approximation domains, respectively fx1 ; . . . ; xN g and fxk ; . . . ; xk‡n g. We denote the radius of the computational domain by R …ˆ xN † and that of the approximation domain G by RG …ˆ xk‡n †. Unless noted otherwise, we assume the following values: risk-free rate of return, r ˆ 0:1 per year; volatility, r2 ˆ …0:3†2 per year; duration, T ˆ 1:0 year; time step size,

M.D. Marcozzi et al. / Appl. Math. Comput. 124 (2001) 197±214

205

l ˆ Ds ˆ 0:045; mesh spacing, h ˆ xi xi 1 ˆ 0:0625. We denote the exact (non-dimensional) put value by p and its numerical approximation by pn . 3.1. Mesh size and arti®cial boundary conditions We consider here the e€ect of space and time step size, computational domain, and arti®cial boundary conditions on the asymptotic performance of the RBF and FD schemes for the European case. For the one-asset case arti®cial boundary conditions are applied by imposing the asymptotic constraints P …S; t† ! Ee P …S; t† ! 0

r…T t†

‡

as S ! 0‡ ;

…3:1a†

as S ! 1

…3:1b†

for all 0 < t < T , once non-dimensionalized, directly onto R and R, respectively. For a uniformly spaced mesh fx1 ; x2 ; . . . ; xN g of size h ˆ xi xi 1 , consecutively ordered, such that x1 ˆ R and xN ˆ R, and a time step of l ˆ Ds, the following one-sided second-order space, implicit ®rst-order time di€erencings were employed at the computational boundaries for the ®nite di€erence approximation: p…x1 ;…m ‡ 1†l†

p…x1 ; ml† ˆ

l ‰2p…x1 ;…m ‡ 1†h† 5p…x2 ;…m ‡ 1†h† h2 ‡ 4p…x3 ; …m ‡ 1†h† p…x4 ;…m ‡ 1†h†Š …3:2a†

and p…xN ; …m ‡ 1†l†

p…xN ; ml† ˆ

l ‰ 12h2

9p…xN ; …m ‡ 1†h† ‡ 7p…xN 1 ; …m ‡ 1†h† ‡ 16p…xN 2 ; …m ‡ 1†h† ‡ 2p…xN 3 ; …m ‡ 1†h†Š

…3:2b†

for m ˆ 1; 2; . . . ; M. 3.2. Domain decomposition and mesh scaling In Figs. 1 and 2, we consider the e€ect of mesh size and computational domain on the RBF and FD approximations. We note, however, for the given parameters, that the RBF approximation e€ectively achieves asymptotic convergence for R=RG  6 whether arti®cial boundary conditions are employed or not; that is, for an approximating domain 2 with RG ˆ 0:2, computationally, 2

An approximation domain of RG ˆ 0:2 constitutes approximately 20% deviations from the strike price.

206

M.D. Marcozzi et al. / Appl. Math. Comput. 124 (2001) 197±214

Fig. 1. E€ect of mesh size h at expiry; RG ˆ 0:1. Note: NBC denotes ``No Boundary Conditions'' while BC indicates that (3.1) has been employed as boundary conditions. (a) RBF (NBC); (b) RBF (BC); (c) FD (NBC); (d) FD (BC).

Fig. 2. E€ect of the computational domain R; RG ˆ 0:2.

the (non-dimensional) boundary ``at in®nity'' may be placed at R ˆ 1:25. This contrasts sharply with results obtained by ®nite-di€erencing where, for a ®xed ratio R=RG , numerical stability appears to depend on a critical mesh size; presumably, inaccuracies due to end e€ects manifest themselves as numerical instabilities with further mesh re®nements. The present arti®cial boundary

M.D. Marcozzi et al. / Appl. Math. Comput. 124 (2001) 197±214

207

conditions are seen to always limit the accuracy of the RBF computations. When employed for R=RG P 6, however, the FD scheme takes on the accuracy and stability of RBF method. In Fig. 3 we note the expected impact of time step size on the RBF and FD calculations and the initial transient. In particular, numerical error present due to the irregularity of the expiry payo€ initially is propagated through the computational domain until ®nally overcome by the compounding e€ects of temporal approximation. Employing methods developed in [12], we may choose the ``optimal'' time step l ˆ Ds such that jpm‡1 pm j ˆ const which we evaluated at the exercise price. We remark that the localized computational domain and adaptive time scaling are key advantages that variational method evidence over lattice based techniques. For valuations encompassing a broad enough expanse, so-called global valuations, the asymptotic behavior (3.1b) will be observed in the put price. This information may be exploited to minimize the e€ort associated with such computations through mesh scaling; the prescription of nodal density based upon the far-®eld behavior of the solution (cf. [15]). For FD approximations, mesh scaling necessarily involves domain decomposition; the coupled resolution of the deterministic put valuations (2.1a), (2.1b) and (2.7a), (2.7b) onto

Fig. 3. E€ect of time step size; R ˆ 1:5, RG ˆ 0:2. (a) Radial Basis Functions; (b) ®nite di€erencing; (c) adaptive error control …h ˆ 0:0097†; (d) adaptive time step size.

208

M.D. Marcozzi et al. / Appl. Math. Comput. 124 (2001) 197±214

exhausting subdomains. For RBF computations, domain decomposition becomes necessary when ill-conditioning e€ects the accuracy of the computation. This may occur, for example, when mesh spacings di€er by orders of magnitude. Domain decomposition algorithms may generally be classi®ed as either overlapping or non-overlapping (cf. [9,11]). We consider here overlapping domains such that the put's price was matched at the interface and the adjacent node. That is, relative to Fig. 4(a), the put value on subdomain 1 at node x5 was equated to node y1 on subdomain 2, similarly at the pairing …x6 ; y2 †. For a uniform mesh, such a construction will have no e€ect on the FD scheme. For the RBF computations, pairings of additional points did not e€ect the accuracy of the calculations to any noticeable extent within the range of parameters considered and so are not presented here (cf. [24]). Mesh scaling with a scale factor of r is illustrated schematically in Fig. 4(b), where for the FD scheme the value of x6 ˆ r 1 y1 ‡ …1 r 1 †y2 . A scale factor of r ˆ 2 was employed throughout the quoted results and the shape parameter e c was calculated for each subdomain. The impact of domain decomposition and mesh scaling can be seen in Fig. 5. In particular, the RBF method remains una€ected while the lower-order coupling of the subdomains for the FD computations appears to ultimately e€ect the accuracy of the method. Fig. 6 represents the global European and American valuations when domain decomposition and mesh scaling are employed with the RBF method. The FD calculations could not be discerned within the resolution of the ®gure and so were not included. Finally, in Fig. 7, we consider the error at expiry associated with the socalled Bermuda approximation to the American formulation. In particular, for a ®xed spatial resolution, we see that as the time step decreases the Bermuda approximation (denoted by pB ) asymptotically approaches the American linear complementary formulation given by (2.11a), (2.11b) and (2.12a), (2.12b).

Fig. 4. Domain decomposition schematic. (a) RBF overlapping subdomains; (b) FD non-overlapping subdomains.

M.D. Marcozzi et al. / Appl. Math. Comput. 124 (2001) 197±214

Fig. 5. E€ect of domain decomposition and mesh scaling.

Fig. 6. European and American put valuations.

Fig. 7. Bermuda approximation.

209

210

M.D. Marcozzi et al. / Appl. Math. Comput. 124 (2001) 197±214

4. The two-asset case Analogous to the analysis of the preceding sections, we consider a (nontrivial) Black±Scholes economy with a ®xed risk-free rate of return r such that the asset prices S1 and S2 satisfy the stochastic system dS1 ˆ rS1 dt ‡ r1 S1 dB1 ; dS2 ˆ rS2 dt ‡ r2 S2 dB2 ; where t represents time, r2i the volatility of asset Si , and dBi Brownian motion with respect to the risk-neutral measure with q ˆ corr …S1 ; S2 †. Supposing a payo€ h such that h…t† ˆ maxfmin…S1 ; S2 †; 0g, the arbitrage-free price of the put in the economy at time t is then  P …S1 ; S2 ; t† ˆ inf E e r…T t† h…t† j Ft ; 06t6T

where again fFt g0 6 t 6 T denotes the ®ltration associated with the Brownian motion and E the expectation with respect to the risk-neutral measure. Formally, the value function satis®es the variational inequality:  oP 1 2 2 o2 P o2 P 1 2 2 o2 P r1 S1 2 qr1 r2 S1 S2 rS min ot 2 oS1 oS2 2 1 1 oS22 oS1 …4:1a†  oP oP rS1 rS2 ‡ rP ; P maxfmin…S1 ; S2 †; 0g oS1 oS2 for Si > 0 and t 2 …0; T †, such that P …S1 ; S2 ; T † ˆ maxfmin…S1 ; S2 †; 0g

…4:1b†

for Si > 0. We may obtain the canonical form of (4.1a) and (4.1b) by making the following substitutions: Let ! 1 r1 S1 ˆ K exp p x2 ; 1 q2 r2 S2 ˆ K exp tˆT

q p x2 1 q2

s …1=2†r22 …1

q2 †

! x1 ;

;

P ˆ K exp …a1 x1 ‡ a2 x2 ‡ bs†p…x1 ; x2 ; s†

M.D. Marcozzi et al. / Appl. Math. Comput. 124 (2001) 197±214

such that a1 ˆ a=2;

a2 ˆ b=2;

where aˆ

……1=2†r21

bˆ

1 2 …a ‡ b2 † 4

r†q…r2 =r1 † ……1=2†r22 …1=2†r22 …1 q2 †

r†

211

c;

;

r† r2 p2 1 q; q2 † r1

……1=2†r21 …1=2†r22 …1 r cˆ 2 …1=2†r2 …1

bˆ

q2 †

in which case we have ( op o2 p o2 p ; p exp ‰ …a1 x1 ‡ a2 x2 ‡ bs†Š min os ox21 ox22 " ! 1 r1 q  min exp p x2 ; exp p x2 1 q2 r2 1 q2

!!#‡ ) x1

ˆ0 …4:2a†

2

for all …x1 ; x2 ; s† 2 R  R‡ , such that " ! 1 r1 p…x1 ; x2 ; 0† ˆ min exp p x2 ; exp 1 q2 r2

q p x2 1 q2

!!#‡ x1 …4:2b†

2

for all …x1 ; x2 † 2 R . We may then de®ne the linear complimentary form and discretization associated with (4.2a) and (4.2b) in the same manner as in

Fig. 8. Two-asset convergence.

212

M.D. Marcozzi et al. / Appl. Math. Comput. 124 (2001) 197±214

Fig. 9. Two-asset value function. (a) Expiry; (b) time t ˆ 0.

Section 2.2. We consider now the approximation of a put option on the minimum of two-assets (cf. [17,26]). Fig. 8 displays the e€ect of mesh size and computational domain on the RBF approximation when no boundary conditions are employed. We assume here that r ˆ 0:10, r1 ˆ r2 ˆ 0:3, q ˆ corr …S1 ; S2 † ˆ 0:0, and RG ˆ 0:2 and note evidently the impact of the computational domain on the convergence of the method. Fig. 9 displays the European put value function.

5. Conclusions The deterministic variational formulations of mathematical ®nance typically involve large numbers of underlying assets (i.e., independent exogenous sources of uncertainty) and the derivation of arti®cial boundary conditions associated with these formulations remain for many problems a substantial, if not insurmountable, diculty (e.g. pricing treasury bonds under an HJM term structure, corporate bonds, options on a foreign currency [1,3]). The results of this paper indicate that such arti®cial boundary conditions are not, in fact, necessary. Moreover, we demonstrate the limitations in accuracies that result from computational domains of insucient radius. If pricing information may be restricted to limited regions of the state space (e.g. a neighborhood of the early exercise boundary), we have shown for two prototype valuations that ``small'' computational domains should suce. These results open the possibility of computations involving a ``large'' number of state variables. To complement this analysis, RBF approximations hold out the promise of more ecient approximations in higher dimensions than tra-

M.D. Marcozzi et al. / Appl. Math. Comput. 124 (2001) 197±214

213

ditional PDE based techniques such as ®nite di€erencing and ®nite elements. Coupled with advances in time step management, it appears that an ecient deterministic variational approach to pricing more general securities may be attainable. References [1] K.I. Amin, R.A. Jarrow, Pricing foreign currency options under stochastic interest rates, J. Int. Money Finance 10 (1991) 310±329. [2] K. Amin, R. Jarrow, Pricing options on risky assets in a stochastic interest rate economy, in: L. Hughston (Ed.), Vasicek and Beyond, Approaches to Building and Applying Interest Rate Models, RISK Pub., 1996, pp. 235±252. [3] K.I. Amin, J.N. Bordurtha Jr., Discrete-time valuation of American options with stochastic interest rates, Rev. Financial Studies 8 (1995) 193±234. [4] A. Bensoussan, On the theory of option pricing, Acta Appl. Math. 2 (1984) 139±158. [5] A. Bensoussan, J.L. Lions, Applications of Variational Inequalities in Stochastic Control, North-Holland, Amsterdam, 1982. [6] P.P. Boyle, J. Evnine, S. Gibbs, Numerical evaluation of multivariate contingent claims, Rev. Financial Studies 2 (1989) 241±250. [7] M. Broadie, P. Glasserman, A stochastic mesh method for pricing high-dimensional American options, Working Paper, 1997. [8] R.E. Carlson, T.A. Foley, The parameter R2 in multiquadric interpolation, Comp. Math. Appl. 21 (1991) 29±42. [9] T.F. Chan, T.P. Mathew, Domain decomposition algorithms, Acta Num. (1994) 61±143. [10] R. Cottle, J-S. Pang, R.E. Stone, The Linear Complementary Problem, Academic Press, New York, 1992. [11] M.R. Dubal, Domain decomposition and local re®nement for multiquadric approximations. I: Second-order equations in one-dimension, J. Appl. Sci. Comp. 1 (1994) 146±171. [12] K. Eriksson, C. Johnson, Adaptive ®nite element methods for parabolic problems. I: A linear model problem, SIAM J. Numer. Anal. 28 (1991) 43±77. [13] C. Franke, R. Schaback, Convergence order estimates of meshless collocation methods using radial basis functions, Adv. Comp. Math. 8 (1998) 381±399. [14] P.R. Garabedian, Partial Di€erential Equations, Chelsea, New York, 1986. [15] C.I. Goldstein, The ®nite element method with nonuniform mesh sizes for unbounded domains, Math. Comp. 36 (1981) 387±404. [16] P. Jaillet, D. Lamberton, B. Lapeyre, Variational inequalities and the pricing of American options, Acta Appl. Math. 21 (1990) 263±289. [17] H. Johnson, Options on the maximum or the minimum of several assets, J. Fin. Quan. Anal. 22 (1987) 277±283. [18] E.J. Kansa, Multiquadratics ± a scattered approximation scheme with applications to computational ¯uid-dynamics ± I: Surface approximation and partial derivative estimates, Comp. Math. Appl. 19 (1990) 127±145. [19] B. Kamrad, P. Ritchken, Multinomial approximating models for options with k state variables, Management Sci. 37 (1991) 1640±1653. [20] I. Karatzas, On the pricing of American options, Appl. Math. Optim. 17 (1988) 37±60. [21] D. Lamberton, B. Lapeyre, Introduction to Stochastic Calculus Applied to Finance, Chapman & Hall, London, 1996. [22] C.E. Lemke, Bimatrix equilibrium points and mathematical programming, Management Sci. 11 (1965) 681±689.

214

M.D. Marcozzi et al. / Appl. Math. Comput. 124 (2001) 197±214

[23] M.D. Marcozzi, On the approximation of optimal stopping problems with application to ®nancial mathematics, SIAM J. Numer. Anal. 22 (5) (2001) 1865±1884. [24] M.D. Marcozzi, S.C. Choi, C.S. Chen, in: C.S. Chen, C. Brebbia, D. Pepper (Eds.), RBF and Optimal Stopping Problems; An Application to the Pricing of Vanilla Options on One Risky Asset, BETECH99, 1999. [25] A. Ravindran, A computer routine for quadratic and linear programming problems, Comm. Assoc. Comp. Mach. 15 (1972) 818±820. [26] R.M. Stultz, Options on the minimum or the maximum of two risky assets, analysis and applications, J. Fin. Econ. 10 (1982) 161±185. [27] P. Wilmott, S. Howison, J. Dewynne, The Mathematics of Financial Derivatives, Cambridge University Press, Cambridge, 1995. [28] R. Zvan, P.A. Forsyth, K.R. Vetzal, Robust numerical methods for PDE models of Asian options, J. Comp. Fin. 1 (1998) 39±78.