Improved radial basis function methods for multi-dimensional option pricing

In this paper, we have derived a radial basis function (RBF) based method for the pricing of financial contracts by solving the Black–Scholes partial differential equation. As an example of a financial contract that can be priced with this method we have chosen the multi-dimensional European basket call option. We have shown numerically that our scheme is second-order accurate in time and spectrally accurate in space for constant shape parameter. For other non-optimal choices of shape parameter values, the resulting convergence rate is algebraic. We propose an adapted node point placement that improves the accuracy compared with a uniform distribution. Compared with an adaptive finite difference method, the RBF method is 20–40 times faster in one and two space dimensions and has approximately the same memory requirements.     

Error saturation in Gaussian radial basis functions on a finite interval

Radial basis function (RBF) interpolation is a “meshless” strategy with great promise for adaptive approximation. One restriction is “error saturation” which occurs for many types of RBFs including Gaussian RBFs of the form ?(x;α,h)=exp(−α²(x/h)²): in the limit h→0 for fixed α, the error does not converge to zero, but rather to E_S(α). Previous studies have theoretically determined the saturation error for Gaussian RBF on an infinite, uniform interval and for the same with a single point omitted. (The gap enormously increases E_S(α).) We show experimentally that the saturation error on the unit interval, x∈[−1,1], is about 0.06exp(−0.47/α²)‖f‖_∞ — huge compared to the O(2π/α²)exp(−π²/[4α²]) saturation error for a grid with one point omitted. We show that the reason the saturation is so large on a finite interval is that it is equivalent to an infinite grid which is uniform except for a gap of many points. The saturation error can be avoided by choosing α?1, the “flat limit”, but the condition number of the interpolation matrix explodes as O(exp(π²/[4α²])). The best strategy is to choose the largest α which yields an acceptably small saturation error: If the user chooses an error tolerance δ, then .     

Algorithms for cardinal interpolation using box splines and radial basis functions

Summary We describe an algorithm for (bivariate) cardinal interpolation which can be applied to translates of basis functions which include box splines or radial basis functions. The algorithm is based on a representation of the Fourier transform of the fundamental interpolant, hence Fast Fourier Transform methods are available. In numerical tests the 4-directional box spline (transformed to the characteristical submodule of ²), the thin plate spline, and the multiquadric case give comparably equal and good results.     

Preconditioning of elliptic problems by approximation in the transform domain

Preconditioned conjugate gradient method is applied for solving linear systemsAx=b where the matrixA is the discretization matrix of second-order elliptic operators. In this paper, we consider the construction of the trnasform based preconditioner from the viewpoint of image compression. Given a smooth image, a major portion of the energy is concentrated in the low frequency regions after image transformation. We can view the matrixA as an image and construct the transform based preconditioner by using the low frequency components of the transformed matrix. It is our hope that the smooth coefficients of the given elliptic operator can be approximated well by the low-rank matrix. Numerical results are reported to show the effectiveness of the preconditioning strategy. Some theoretical results about the properties of our proposed preconditioners and the condition number of the preconditioned matrices are discussed.     

The truncated Stieltjes moment problem solved by using kernel density functions

In this work the problem of the approximate numerical determination of a semi-infinite supported, continuous probability density function (pdf) from a finite number of its moments is addressed. The target space is carefully defined and an approximation theorem is proved, establishing that the set of all convex superpositions of appropriate Kernel Density Functions (KDFs) is dense in this space. A solution algorithm is provided, based on the established approximate representation of the target pdf and the exploitation of some theoretical results concerning moment sequence asymptotics. The solution algorithm also permits us to recover the tail behavior of the target pdf and incorporate this information in our solution. A parsimonious formulation of the proposed solution procedure, based on a novel sequentially adaptive scheme is developed, enabling a very efficient moment data inversion. The whole methodology is fully illustrated by numerical examples.     

Implementation of a modified Marder–Weitzner method for solving nonlinear eigenvalue problems

Our goal is to propose four versions of modified Marder–Weitzner methods and to present the implementation of the new-type methods with incremental unknowns for solving nonlinear eigenvalue problems. By combining with compact schemes and modified Marder–Weitzner methods, six schemes well suited for the calculation of unstable solutions are obtained. We illustrate the efficiency of the new algorithms by using numerical computations and by comparing them with existing methods for some two-dimensional problems.     

Preconditioners for spectral element methods for elliptic and parabolic problems

In this paper we propose preconditioners for spectral element methods for elliptic and parabolic problems. These preconditioners are constructed using separation of variables and are easy to invert. Moreover they are spectrally equivalent to the quadratic forms which they are used to approximate.     

Convergence and optimization of the parallel method of simultaneous directions for the solution of elliptic problems

For the solution of elliptic problems, fractional step methods and in particular alternating directions (ADI) methods are iterative methods where fractional steps are sequential. Therefore, they only accept parallelization at low level. In [T. Lu, P. Neittaanmäki, X.C. Tai, A parallel splitting-up method for partial differential equations and its applications to Navier–Stokes equations, RAIRO Modél. Math. Anal. Numér. 26 (6) (1992) 673–708], Lu et al. proposed a method where the fractional steps can be performed in parallel. We can thus speak of parallel fractional step (PFS) methods and, in particular, simultaneous directions (SDI) methods. In this paper, we perform a detailed analysis of the convergence and optimization of PFS and SDI methods, complementing what was done in [T. Lu, P. Neittaanmäki, X.C. Tai, A parallel splitting-up method for partial differential equations and its applications to Navier–Stokes equations, RAIRO Modél. Math. Anal. Numér. 26 (6) (1992) 673–708]. We describe the behavior of the method and we specify the good choice of the parameters. We also study the efficiency of the parallelization. Some 2D, 3D and high-dimensional tests confirm our results.     

Some functional equations in the spaces of generalized functions

Summary. Making use of the fundamental solution of the heat equation we find the solutions of some functional equations such as the Cauchy equations, Pexider equations, quadratic functional equations and dAlembert equations in the spaces of Schwartz distributions and Sato hyperfunctions.     

<Emphasis Type="Italic">L</Emphasis><Superscript>1</Superscript>-minimization methods for Hamilton–Jacobi equations: the one-dimensional case

A new approximation technique based on L ¹-minimization is introduced. It is proven that the approximate solution converges to the viscosity solution in the case of one-dimensional stationary Hamilton–Jacobi equation with convex Hamiltonian. This material is based upon work supported by the National Science Foundation grant DMS-0510650. J.-L. Guermond is on leave from LIMSI, UPRR 3251 CNRS, BP 133, 91403 Orsay Cedex, France.     

The fast Fourier transform and the numerical solution of one-dimensional boundary integral equations

Summary Here we present a fully discretized projection method with Fourier series which is based on a modification of the fast Fourier transform. The method is applied to systems of integro-differential equations with the Cauchy kernel, boundary integral equations from the boundary element method and, more generally, to certain elliptic pseudodifferential equations on closed smooth curves. We use Gaussian quadratures on families of equidistant partitions combined with the fast Fourier transform. This yields an extremely accurate and fast numerical scheme. We present complete asymptotic error estimates including the quadrature errors. These are quasioptimal and of exponential order for analytic data. Numerical experiments for a scattering problem, the clamped plate and plane estatostatics confirm the theoretical convergence rates and show high accuracy.     

Fourier analysis of schwarz alternating methods for piecewise hermite bicubic orthogonal spline collocation

The rates of convergence of two Schwarz alternating methods are analyzed for the iterative solution of a discrete problem which arises when orthogonal spline collocation with piecewise Hermite bicubics is applied to the Dirichlet problem for Poisson's equation on a rectangle. In the first method, the rectangle is divided into two overlapping subrectangles, while three overlapping subrectangles are used in the second method. Fourier analysis is used to obtain explicit formulas for the convergence factors by which theH ¹-norm of the errors is reduced in one iteration of the Schwarz methods. It is shown numerically that while these factors depend on the size of overlap, they are independent of the partition stepsize. Results of numerical experiments are presented which confirm the established rates of convergence of the Schwarz methods.This research was supported in part by funds from the National Science Foundation grant CCR-9103451.     

A posteriori error estimators for the Stokes equations II non-conforming discretizations

Summary We present an a posteriori error estimator for the non-conforming Crouzeix-Raviart discretization of the Stokes equations which is based on the local evaluation of residuals with respect to the strong form of the differential equation. The error estimator yields global upper and local lower bounds for the error of the finite element solution. It can easily be generalized to the stationary, incompressible Navier-Stokes equations and to other non-conforming finite element methods. Numerical examples show the efficiency of the proposed error estimator.     

Discrete mechanics and its application to the solution of the n-body problem

A theory of discrete mechanics is developed based on the results of D. Greenspan. Discrete dynamical equations in an inertial frame, in a coordinate system related to some material point, and in a rotating frame are given and the consistency, stability, and convergence of the methods are studied and some numerical examples presented.     

Circulant preconditioners for pricing options

We use the normalized preconditioned conjugate gradient method with Strang’s circulant preconditioner to solve a nonsymmetric Toeplitz system A_nx=b, which arises from the discretization of a partial integro-differential equation in option pricing. By using the definition of family of generating functions introduced in [16], we prove that Strang’s circulant preconditioner leads to a superlinear convergence rate under certain conditions. Numerical results exemplify our theoretical analysis.     

Direct and iterative solution of the generalized Dirichlet–Neumann map for elliptic PDEs on square domains

In this work we derive the structural properties of the Collocation coefficient matrix associated with the Dirichlet–Neumann map for Laplace’s equation on a square domain. The analysis is independent of the choice of basis functions and includes the case involving the same type of boundary conditions on all sides, as well as the case where different boundary conditions are used on each side of the square domain. Taking advantage of said properties, we present efficient implementations of direct factorization and iterative methods, including classical SOR-type and Krylov subspace (Bi-CGSTAB and GMRES) methods appropriately preconditioned, for both Sine and Chebyshev basis functions. Numerical experimentation, to verify our results, is also included.     

The generalized Dirichlet–Neumann map for linear elliptic PDEs and its numerical implementation

A new approach for analyzing boundary value problems for linear and for integrable nonlinear PDEs was introduced in Fokas [A unified transform method for solving linear and certain nonlinear PDEs, Proc. Roy. Soc. London Ser. A 53 (1997) 1411–1443]. For linear elliptic PDEs, an important aspect of this approach is the characterization of a generalized Dirichlet to Neumann map: given the derivative of the solution along a direction of an arbitrary angle to the boundary, the derivative of the solution perpendicularly to this direction is computed without solving on the interior of the domain. This is based on the analysis of the so-called global relation, an equation which couples known and unknown components of the derivative on the boundary and which is valid for all values of a complex parameter k. A collocation-type numerical method for solving the global relation for the Laplace equation in an arbitrary bounded convex polygon was introduced in Fulton et al. [An analytical method for linear elliptic PDEs and its numerical implementation, J. Comput. Appl. Math. 167 (2004) 465–483]. Here, by choosing a different set of the “collocation points” (values for k), we present a significant improvement of the results in Fulton et al. [An analytical method for linear elliptic PDEs and its numerical implementation, J. Comput. Appl. Math. 167 (2004) 465–483]. The new collocation points lead to well-conditioned collocation methods. Their combination with sine basis functions leads to a collocation matrix whose diagonal blocks are point diagonal matrices yielding efficient implementation of iterative methods; numerical experimentation suggests quadratic convergence. The choice of Chebyshev basis functions leads to higher order convergence, which for regular polygons appear to be exponential.     

Matrix computations and polynomial root-finding with preprocessing

Solution of homogeneous linear systems of equations is a basic operation of matrix computations. The customary algorithms rely on pivoting, orthogonalization and SVD, but we employ randomized preprocessing instead. This enables us to accelerate the solution dramatically, both in terms of the estimated arithmetic cost and the observed CPU time. The approach is effective in the cases of both general and structured input matrices and we extend it and its computational advantages to the solution of nonhomogeneous linear systems of equations, matrix eigen-solving, the solution of polynomial and secular equations, and approximation of a matrix by a nearby matrix that has a smaller rank or a fixed structure (e.g., of the Toeplitz or Hankel type). Our analysis and extensive experiments show the power of the presented algorithms.     

Static condensation,partial orthogonalization of basis functions,and ILU preconditioning in the hp-FEM

Static condensation of internal degrees of freedom, partial orthogonalization of basis functions, and ILU preconditioning are techniques used to facilitate the solution of discrete problems obtained in the hp-FEM. This paper shows that for symmetric linear (not necessarily positive-definite) problems, under mild technical assumptions, these three techniques are completely equivalent. In fact, the same matrices can be obtained by the same arithmetic operations. The study can be extended to nonsymmetric problems naturally.