Efficient pricing of European options on two underlying assets by frame duality

Abstract

The PROJ method for pricing European options on one underlying asset was proposed by J. Lars Kirkby and then was applied to price Bermudan and Asian options. In this paper, we extend the method to higher dimensions, especially two-dimensions in which some exotic options can be priced. Our method does not rely on a-priori truncation of the integration range and exhibits excellent performance compared with other state-of-the-art methods, particularly for fatter-tailed short maturity models. We also discuss the errors introduced in each approximation and give corresponding error bounds. Numerical results on implementation of this method to price for popular two-assets options, under both the geometric Brownian Motion and Variance-Gamma dynamics, demonstrate remarkable accuracy and robustness.

Introduction

In financial markets, investors use options to speculate or hedge against investment risk. Hence, option pricing, i.e. knowing the fair price of option contract, is one of the major problems both in theoretical research and practical implementations. Ever since the ground-breaking Black-Scholes option pricing model in 1973, a variety of models have been proposed to modify its inadequacy and hence present better fit with observed market features. However, closed-form formulas for option prices are generally unavailable or not easily derived under these advanced models, even for European vanilla options. For this reason, many numerical methods have been proposed to derive approximate solutions of option pricing problems. With the advantage of flexible applicability and simple implementation, the Monte Carlo (MC) simulation method is widely used to price for financial derivatives at the expense of time, and is especially suitable for pricing path dependent options and high-dimensional options. Different approaches are finite difference methods and binomial tree methods, interested readers are referred to [28] [2].

In this paper, we focus on Fourier inversion based methods. By taking the Fourier transform of the exponentially damped option price (or option payoff), Carr and Madan [3] express the option price as a Fourier transformed integral and then apply the fast Fourier transform (FFT) algorithm to speed up the computation. The FFT method is widely adopted and has since then been adapted to a variety of option payoffs and different models as long as the characteristic function is available. The accuracy of FFT depends on the choice of the damping coefficient. Fang and Oosterlee [11] develop an option pricing method for European options based on Fourier-cosine series expansions, which is called the COS method. This method relies on a priori truncation of the density function range prescribed by cumulants and generates solutions that are very sensitive to the choice of truncation range. Once an appropriate integration interval is determined, the COS method can achieve a high accuracy in milliseconds and has wide applicability to a variety of exotic options under different pricing models. Later, the COS method is used to price for early-exercise and discrete barrier options and European-style Asian options under exponential Lévy models (see [12] [30]), and for Bermudan and barrier options under Heston's model in [13].

To remedy the deficiency of the COS method, local wavelets bases have been considered to recover the underlying density function from its characteristic function instead of from the cosine series expansion. Among them, Haar and B-spline wavelets with linear order are considered in [24] with density coefficients computed via the Cauchy's integral theorem, and higher order B-spline bases are introduced to option valuation by Kirkby [16] with density coefficients computed by Parseval's identity. The option pricing method of density projection (onto B-Splines) by frame duality is termed as PROJ method and was applied to price Bermudan and Asian options (see [22] [17]). Because of the simple form of the compactly supported cardinal B-splines, option pricing by the PROJ method presents accuracy and efficiency even for heavy tailed asset processes. The SWIFT method, based on the Shannon wavelet expansion of the underlying density function, has been applied for pricing European options [25]. The SWIFT technique is efficient for the valuation of long- and short-maturity option contracts and it doesn't rely on a priori truncation of computation interval, since it can automatically calculate the number of terms in series expansion to meet a predefined tolerance error. When analytical expression for the payoff coefficients is not available, closed form quadrature rules can be applied, e.g. the Newton-Cotes quadrature formulas suggested in [18] and further illustrated in [19] for barrier options, as well as in [21] with complex payoffs in the case of swing option pricing. It should be noted that even if the exact integral of payoff coefficients exists, it may lead to numerical instability when the basis resolution is refined, as discussed in detail in [18]. In this case, closed form approximations to the coefficient functions can serve as a competitive alternative to produce more stable results, which also enable the use of higher order (and faster converging) B-spline scaling functions, see [18] for more details.

The previous mentioned numerical option pricing methods are applicable for European-style options with one underlying asset driven by a single stochastic process. There is a class of exotic options with multiple underlying assets, called the multicolor rainbow options, whose payoffs depend on the performance of multiple underlying assets at maturity. For valuation of exotic options on two underlying assets, closed-form pricing formulas are unavailable under various models with few exceptions, therefore there is a high demand for efficient numerical methods to give accurate option prices. For instance, a analytical pricing formula exists for geometric basket options under the geometric Brownian motions (GBM). In [26], the COS method is extended to higher dimensions to be able to price options on two or more underlying assets, referred to as the 2D-COS method, which is a highly efficient method particularly for exponential Lévy models. Analogously, the SWIFT method is also generalized to multidimensional case [7] and is called 2D-SWIFT method in two-dimensional case. The advantage of 2D-SWIFT compared with the 2D-COS method is the same as in the one-dimensional case. In multidimensional option pricing, especially in the two-dimensional case, the fast Fourier transform (FFT) approach of Carr and Madan is extended by Hurd and Zhou [15] to be able to calculate the fair price for spread options. Spread options are difficult to price without the use of Monte Carlo simulation or the 2D-COS method. By relating the option price with a Fourier transformed integral through a novel Fourier representation of the basic spread payoff function $C(x_1,x_2)={(e^{x_1}-e^{x_2}-1)}^{+}$, where $x_i=\ln (S_T^i)$ and the strike price $K=1$, they developed a two-dimensional version of the valuation formula of Carr and Madan. In fact, this method can be applied to options whose complex Fourier transform of the payoff exists when imposing restrictions on the imaginary part in frequency domain.

In this paper, we extend the PROJ method to higher dimensions to be capable of pricing multicolor rainbow options and we focus on the two-dimensional case, called the 2D-PROJ method. The density function in the risk-neutral pricing formula is approximated by its projection onto a closed subspace of $L^2$ and formulated as a finite combination of multidimensional B-spline scaling functions. With the help of frame duality theory, the projection coefficients of density function, also referred to as density coefficients, have an exact integral representation. Moreover, the integral of product of the payoff function and B-spline basis function, called the payoff coefficient, either has analytical expressions for some two-color rainbow options or can be calculated numerically. The motivation behind this paper is that we find the analytical payoff coefficients produced by mathematical modeling software are too lengthy to be used in implementation of the 2D-COS method. Furthermore, they are difficult to be calculated by hand. However, benefiting from the simple form and local behavior of cardinal B-splines, we can easily derive an analytical expression for payoff coefficients. In addition to that, we do not rely on a-priori integration interval. For highly peaked densities, our method can generate robust prices within milliseconds if parameters are suitably selected. Although approximations occur several times in the 2D-PROJ algorithm, we discuss the errors produced in each approximate step and confirm an exponential rate of convergence. From a practical point-of-view, the 2D-PROJ method outperforms its nearest rival, 2D-SWIFT [7], in the context of easy to be implemented since there are only two parameters in this approach.

The remainder of this paper is organized as follows. We start in Section 2 by presenting the option pricing and providing a brief review of multiresolution analysis in higher dimensions and B-spline scaling functions. In Section 3 the 2D-PROJ pricing formula is derived and extended to higher dimensions. An error analysis is included in Section 4. We perform a series of numerical results to demonstrate the robustness and accuracy of the 2D-PROJ method in Section 5. The advantage of the 2D-PROJ method is discussed in Section 6. Finally, the paper is concluded in Section 7.