An Algorithm for Linear Complementarity and its Application in American Options Pricing

Abstract

In this talk, we consider the numerical solution of American options pricing problems. Most options traded on options exchanges world-wide and a large fraction of options traded over-the-counter are of the American-style, including options on stocks of individual companies, stock indexes, foreign currencies, interest rates, commodities, and energy. Options books of a large financial institution may contain options on thousands of different underlying assets, and perhaps several dozen different contracts (with expiration dates ranging from days to years, and different strike prices). As the underlying asset prices change throughout the trading day, the options prices change as well. Re-pricing a large options book in real time may thus require re-computing thousands of options prices quickly. For such large scale applications, fast numerical algorithms are essential. When the prices of underlying assets are assumed to follow a diffusion process, such as in the classical BlackScholes-Merton model based on the geometric Brownian motion process, or in extensions such as Heston's stochastic volatihty model, the pricing function of an American-style option satisfies a system of parabolic partial differential variational inequahties. After this system is discretized in space and time, it yields a linear complementarity problem, which must be solved at each time step. Thus, the fast solution of linear complementarity problems (LCPs) is of great practical importance in computational finance. The most popular LCP method at present is the projected SOR iteration, or the closely related variant, the projected Gauss-Seidel iteration [2]. The standard treatment of LCPs for American option pricing can be found, for example, in [8] for the simple case of the Black-Scholes-Merton model and in [4] for several more complicated settings. Several new active-set methods [1,7] have recently been proposed for solving these LCPs more efficiently. Some of the most promising results are reported by Borici and Luethi [1], who developed a variant of the simplex-like method for LCPs with Z-matrices [2]. We argue in this talk that much greater speedups can be obtained with an algorithm that combines iterations of the projected Gauss-Seidel (or SOR) method with reduced-space steps. This two-phase approach exploits the fact that the projected Gauss-Seidel iteration often makes a quick estimation of the optimal active set, while the reduced-space iteration can dramatically improve upon this estimate and yield a fast rate of convergence. We illustrate the performance of this algorithm on both the Black-Scholes-Merton model (using various values of volatility and maturity) and the Heston model [3] with stochastic volatility. Let us begin by describing the Black-Scholes-Merton model. Consider an American put option with strike price K > 0 and maturity time T > 0. If the option is exercised when the underlying asset price is S, the option holder receives the payoff W{S) = {K -S)^ = max{K -S,0). Similarly, the payoff function for an American call option is ^{S) = {SK)+. Let V{t, S) be the option value at time t G [0, T] when the asset price is S. We assume that V solves the following partial differential variational inequahty (see, e.g., [5]):