Sum of all Black-Scholes-Merton models: An efficient pricing method for spread, basket, and Asian options

Introduction

Ever since the celebrated success of the Black-Scholes-Merton (BSM) model, the effort to extend its simple analytic solution to the derivatives written on multiple underlying assets has been of great interest to researchers (Broadie and Detemple 2004). One important class of such derivatives is the options on a linear combination of assets following correlated geometric Brownian motions (GBMs); these include the following three popular option types:

Spread option: European-style option on the difference between two asset prices;

Basket option: European-style option on the sum of multiple asset prices with positive weights;

Asian option: option on the average price of one underlying asset on a pre-determined discrete time set or continuous time range.

These three option types arguably represent the most actively traded non-vanilla options on exchanges or in over-the-counter markets. This is because a linear combination is a common way of associating multiple prices – of different assets or various times – and such options provide customized hedge or risk exposure. For examples and financial motivations, see the introductions of Carmona and Durrleman 2003 and Linetsky 2004.

However, such option pricings under the BSM model, not to mention the models beyond the BSM, is not a trivial matter. This is because, unlike normal random variables (RVs), a linear combination of correlated log-normal RVs neither falls back to the same class of distribution nor has a distribution expressed in any analytic form in general. Therefore, the exact valuation of option prices involves a multidimensional integral over positive payoffs under the risk-neutral measure. The numerical evaluation of such an integral, however, suffers from the curse of dimensionality. For example, even the coarse discretization of a standard normal distribution from $-5$ to $5$ with a grid size of $0.25$ (41 points per dimension) leads to 3 million points for four assets and 116 million points for five assets, substantially exceeding the size of a typical Monte Carlo simulation.

2. Literature review

A vast amount of the literature has examined each of the three option types. First, a review is conducted of the analytic methods of approximating the option price or the lower and upper bounds of the option price. For the spread option, Kirk’s formula (Kirk 1995) that is widely used in practice is an approximate generalization of Margrabe’s formula (Margrabe 1978) for an exchange option, that is, a spread option with zero-strike price; several improvements to the formula have followed (Bjerksund and Stensland 2014; Lo 2015). Carmona and Durrleman 2003 computes the lower bound of a spread option price as the maximum over the prices from all possible linear, and therefore sub-optimal, exercise boundaries. Li et al. 2008 proposes a closed-form formula based on quadratic approximation of the exercise boundary.

Basket and Asian options share the underlying ideas for analytic approximation because the payoff of an Asian option depends on a basket of correlated prices over the observation period. One of the most popular approaches is to approximate the distribution of the GBM sum with other analytically known distributions, such as log-normal (Levy and Turnbull 1992; Levy 1992), reciprocal gamma (Milevsky and Posner 1998a; Milevsky and Posner 1998b), shifted log-normal (Borovkova et al. 2007), and log-extended-skew-normal (Zhou and Wang 2008) distributions, and with perturbation expansions from known distributions (Turnbull and Wakeman 1991; Ju 2002). Another popular idea is to exploit the geometric mean of GBMs—as opposed to arithmetic mean—whose distribution is log-normal and hence analytically solvable. The option price on the geometric mean is a reasonable proxy for that on the arithmetic mean (Gentle 1993). Thus, it can be used as a control variate, thereby reducing the Monte Carlo variance for Asian (Kemna and Vorst 1990) and basket (Krekel et al. 2004) options. Curran 1994 uses the geometric mean as a conditioning variable to analytically estimate option prices. The conditioning approach is further refined by Beisser 1999 and Deelstra et al. 2004 and also applied to the continuously monitored Asian options (Rogers and Shi 1995). For other basket and Asian option pricing approaches as well as classifications, see Zhou and Wang 2008 and the references therein.

Analytic approximation methods are appealing because of their simple computation, but has a major limitation in that the results are not precise. While each method is accurate for certain parameter ranges where the underlying assumptions are valid, the accuracy of any one method can hardly be validated for all ranges. For example, no single method performs well in basket options under various parameter sets (Krekel et al. 2004), although the method of Ju 2002 is outstanding overall. Therefore, practitioners must carefully identify the parameter range in which the method of interest works best. Given the multi-asset aspect of such problems, the charting of parameter maps in advance is not a trivial matter. Moreover, since errors cannot be controlled in approximation methods, the tendency is to eventually apply external methods, typically a Monte Carlo simulation, to obtain a benchmark value.

Convergent pricing methods are fewer in number compared to approximation methods. Convergent implies that the method can produce a deterministic price, as opposed to Monte-Carlo methods, and converge to the true value with a reasonable amount of computation when the computational parameters (e.g., grid size) are tuned. Convergent methods are feasible for spread options because the problem is two-dimensional. Ravindran 1993 and Pearson 1995 reduce the pricing problem to a one-dimensional integration over the BSM prices with varying spot and strike prices. In addition, Dempster and Hong 2002 and Hurd and Zhou 2010 apply a two-dimensional fast Fourier transform (FFT).

To the best of the author’s knowledge, few studies on basket options use methods that can be defined as convergent. Particularly, no attempt has been made to use direct integration, even for error measurement, thus indicating the challenge posed by the approach (see § 4.1). Although the FFT approach proposed by Leentvaar and Oosterlee 2008 reduces the computation time through parallel partitioning, it does not significantly reduce the computation amount.

Previous convergent methods used for Asian options are based on the idea that pricing involves a single price process over time. The continuously averaged Asian option, although hardly traded in practice owing to contractual difficulty, has analytic solutions comprising the triple integral (Yor 1992), Laplace transform in maturity (Geman and Yor 1993), and a series expansion (Linetsky 2004). For discrete averaging, a series of studies have exploited the recursive convolution, referred to as the Carverhill-Clewlow-Hodges factorization, on the probability density function (PDF) or price (Carverhill and Clewlow 1990; Benhamou 2002; Fusai and Meucci 2008; Černỳ and Kyriakou 2011; Fusai et al. 2011; Zhang and Oosterlee 2013). Cai et al. 2013 describes the price of a discretely monitored Asian option as an asymptotic expansion on a small observation interval.

Despite structural similarity, only a few studies are applicable to all the three option types. Carmona and Durrleman 2005 extend the lower bound approach (Carmona and Durrleman 2003) to a multi-asset problem, while Deelstra et al. 2010 use the commonality theory to approximate the option price. However, because these studies follow the analytic approximation class, their methods have the aforementioned limitations. Essentially, no convergent method can work consistently for all the three option types.

3. Contribution of paper

The aim of this study is to find an efficient and unified pricing method for options on the linear combinations of correlated GBMs. While the common assumption is that outright multidimensional integration is computationally prohibitive, an innovative integration scheme is developed, which can significantly alleviate the curse of dimensionality. First, it is observed that if all the price processes are driven by a single Brownian motion, the option price can be analytically obtained by a multidimensional extension of the BSM formula, with the exercise boundary obtained from a numerical root-finding. Therefore, the option price can be integrated analytically for the first dimension and numerically for the remaining dimensions. The key to the current approach is to choose the first dimension via factor rotation such that the analytic price from the first dimension behaves well (i.e., smooth and slowly varying) as an integrand for the numerical integration that follows. Numerical experiments show that even a coarse discretization produces a very accurate price, and that the price quickly converges to the true value as the number of nodes increase.

The contributions of this study in the context of each option type are presented in order below. Here, the approach to spread options is similar to those of Ravindran 1993 and Pearson 1995; that is, the integration in the first dimension is done analytically. However, the factor rotation in this study’s method significantly reduces the cost of numerical integration in the second dimension. Even at the expense of extra computation due to numerical root-finding, which is not found in the aforementioned studies, the overall computation is much lighter owing to reduced discretization on the second dimension.

As regards basket options, the method used here is fully convergent. While the speed of convergence depends on the covariance structure, in general, the method can converge to the true option value for a wide range of parameters and dimensions. For the first time, the converged prices for several benchmark tests used in the literature can be reported.

As for the study on Asian options, this study’s pricing scheme serves as a novel alternative to those of previous works. The integration approach is particularly problematic for Asian options owing to large dimensionality, that is, a large number of observations. In effect, the study utilizes several of the first factors of a Brownian motion series representation as with the principal component analysis (PCA). Therefore, the method falls short of being truly convergent. However, the method is accurate for all practical purposes and enables cheaper computation compared to existing methods. Furthermore, it is capable of pricing a continuously monitored Asian option in a discrete monitoring framework, which is a rare approach compared to those in the opposite direction. To this extent, the method can be used to flexibly handle features such as non-uniform weights, non-uniform averaging intervals (e.g., forward-start averaging), and time-dependent volatility, which are difficult to incorporate in methods based on the continuum theory (Linetsky 2004; Cai et al. 2013; Fusai et al. 2011). Asian options are discussed in detail in § 5.

This study focuses on the BSM model, but its findings can be useful for other models too. The results can be trivially modified to apply to displaced GBMs, as shown in a numerical example in § 6. Displaced GBMs have an extra degree of freedom to capture volatility skew, if not full smile, as observed in the option market. The results can also be applied to stochastic volatility models, such as the Hull and White 1987, Heston 1993, and stochastic-alpha-beta-rho (Hagan et al. 2002) models. In these models, the options are priced using the conditional Monte-Carlo method, where the price is expressed as an expectation of the BSM prices, conditional on quantities such as terminal volatility and integrated variance; see Willard 1997, Broadie and Kaya 2006, and Cai et al. 2017, respectively. Thus, an efficient BSM pricing method is critical for pricing the spread or basket options under such approaches. However, a detailed implementation is beyond the scope of this study.

This paper is organized as follows. Section 2 formulates the problem. Section 3 outlines the multidimensional integration scheme. Section 4 discusses the optimal rotation of the dimensions facilitating integration. Section 5 discusses the implications for Asian options. Section 6 reports the numerical results. Finally, Section 7 concludes the study.

Model setup and preliminaries

Assume that asset prices, $S_{k}$ for $1\leq k\leq N$, follow the correlated GBMs under risk-neutral measure

where $\sigma_{k}$ is the volatility, $q_{k}$ is the dividend rate, $r$ is the risk-free interest rate, and $W_{k}(t)$ is a standard Brownian motion with correlation $dW_{k}(t)dW_{j}(t)=\rho_{kj}dt$ ($\rho_{kk}=1$). The final payoff of the options considered here depends on a linear combination of the asset prices observed earlier than or at the expiry $T$. The payoff of a vanilla call option with strike price $K$ is as follows:

for weights, $w_{k}$, and observation times, $0\leq t_{k}\leq T$. Here, $(x)^{+}=\max(x,0)$ is the positive-part operator. The current setup is generic enough to include, but not limited to, the following three option types:

European spread option: $w_{k}<0$ for some, but not all, $k$; additionally, $t_{k}=T$ for all $k$.

European basket option: $w_{k}>0$ and $t_{k}=T$ for all $k$.

Asian option: $w_{k}>0$ for all $k$ with $\sum w_{k}=1$ and $0\leq t_{1}<\cdots<t_{N}=T$. The price processes are all identical, $S_{k}(t)=S_{j}(t)$ for $k\neq j$; thus, the index $k$ is omitted without ambiguity; for example, $S(t)$, $W(t)$, $\sigma$, and $q$. The continuously monitored Asian option, whose payoff is $(\frac{1}{T}\int_{0}^{T}S(t)dt-K)^{+}$, will be considered under the discrete framework.

Next, a few notations and conventions are set. For matrix $\boldsymbol{A}$, the following is noted: the $k$-th row vector of $\boldsymbol{A}$, the $j$-th column vector of $\boldsymbol{A}$, and the $(k,j)$ component of $\boldsymbol{A}$ by $\boldsymbol{A}_{k*}$, $\boldsymbol{A}_{*j}$, and ${A}_{kj}$, respectively. Using these notations, for example, the matrix multiplication $\boldsymbol{C}=\boldsymbol{A}\boldsymbol{B}$ can be expressed as ${C}_{kj}=\boldsymbol{A}_{k*}\boldsymbol{B}_{*j}$. The transpose of $\boldsymbol{A}$ is noted by $\boldsymbol{A}^{\scriptscriptstyle T}$ and the identity matrix by $\boldsymbol{I}$. For vector $\boldsymbol{x}$, the $k$-th component is noted by $x_{k}$. Moreover, unless otherwise stated, the vector is a column vector. The $L^{2}$-norm of $\boldsymbol{x}$ is $|\boldsymbol{x}|=\sqrt{\boldsymbol{x}^{\scriptscriptstyle T}\boldsymbol{x}}$ and the Frobenius norm of $\boldsymbol{A}$ is $\|\boldsymbol{A}\|_{F}=(\sum_{k,j}{A}_{kj}^{2})^{1/2}$. Further, unless otherwise specified, the dimensions of a matrix are $N\times N$, the size of a vector is $N$, and the indices, $k$ and $j$, run from $1$ to $N$. For the factor matrix $\boldsymbol{V}$ to be defined below, $k$ is used for indexing rows (assets) and $j$ for indexing columns (Brownian motions or factors). Throughout this study, the terms factor and dimension are used interchangeably.

Under the BSM model, the log prices, $\log S_{k}(t_{k})$, follow correlated normal distributions, with the covariance matrix $\boldsymbol{\Sigma}$ given as

If $\boldsymbol{V}$ is a square root matrix of $\boldsymbol{\Sigma}$, satisfying $\boldsymbol{V}\boldsymbol{V}^{\scriptscriptstyle T}=\boldsymbol{\Sigma}$, then the observation $S_{k}(t_{k})$ can be decorrelated to

where $\boldsymbol{z}$ is a vector of independent standard normal RVs and $F_{k}=S_{k}(0)\,e^{(r-q_{k})t_{k}}$ is the $t_{k}$-forward price observed at $t=0$. The symbol $\,{\mathrel{\mathop{\kern 0.0pt=}\limits^{d}}}\,$ stands for equality in distribution law. Since $\boldsymbol{V}$ is multiplied to the state vector $\boldsymbol{z}$, it is referred to as a risk factor matrix, or simply a factor matrix. The square root matrix $\boldsymbol{V}$ is not unique. Although the Cholesky decomposition $\boldsymbol{C}$ is a popular choice, any matrix rotated from $\boldsymbol{C}$, such as $\boldsymbol{V}=\boldsymbol{C}\boldsymbol{Q}$ for an orthonormal matrix $\boldsymbol{Q}$, is also a square root matrix of $\boldsymbol{\Sigma}$. However, note that the norm of row vectors $|\boldsymbol{V}_{k*}|$ is invariant under any rotation; since $|\boldsymbol{V}_{k*}|^{2}$ is the variance of the $k$-th asset’s return, $|\boldsymbol{V}_{k*}|^{2}=\boldsymbol{V}_{k*}\boldsymbol{V}_{*k}={\Sigma}_{kk}$. Further, the Frobenius norm of $\boldsymbol{V}$ is also invariant because $\|\boldsymbol{V}\|_{F}^{2}=\sum_{k}|\boldsymbol{V}_{k*}|^{2}=\sum_{k}{\Sigma}_{kk}$.

The forward value of the call option price becomes an $N$-dimensional integration,

where $n(\boldsymbol{z})$ is a multivariate standard normal PDF.

Integration scheme

The BSM model’s analytic tractability can be extended to the multi-asset case when a single Brownian motion drives all assets. Consider the integration over the first dimension $z_{1}$ only:

where the dependence on other dimensions is absorbed into the coefficient function defined as

While $\dot{\boldsymbol{z}}$ has $z_{1}=0$ for computational convenience, $F(\dot{\boldsymbol{z}})$ for function $F$ should be understood as $F(z_{2},\cdots,z_{N})$. Note that $f_{k}(\dot{\boldsymbol{z}})>0$ and $E(f_{k}(\dot{\boldsymbol{z}}))=1$ for all $k$. The value $C_{\text{BS}}(\dot{\boldsymbol{z}})$ can be seen as the price of an option on the weighted asset prices driven by a single Brownian motion, where the forward price is $F_{k}\,f_{k}(\dot{\boldsymbol{z}})$ and the standard deviation of the log price is ${V}_{k1}$ for the $k$-th asset.

This one-dimensional integration can be analytically evaluated only if the region of the positive payoff is identified. However, to find the roots of the payoff, a numerical method must be used, like the Newton-Raphson method. As explained below in § 4.2, $\boldsymbol{V}$ can be chosen such that the payoff monotonically increases in $z_{1}$ and there always exists a unique root, $z_{1}=-d(\dot{\boldsymbol{z}})$, of the equation

The integration from $-d(\dot{\boldsymbol{z}})$ to $\infty$ yields

This is a multi-asset extension of the BSM formula. The original BSM formula is a special case of the single asset case ($N=1$, $t_{1}=T$, and $w_{1}=1$):

where $\boldsymbol{V}$ is a scalar value, ${V}_{11}=\sqrt{{\Sigma}_{11}}=\sigma_{1}\sqrt{T}$.

Despite the cumbersome numerical root-finding, the analytic integration of $z_{1}$ plays an important role, which involves more than simply reducing one dimension in the integration. First, because of the cusp of the option payoff at the strike, numerical integration on the first factor would have required the densest discretization. Thus, analytic integration allows us to skip the most computationally costly dimension, albeit at the expense of numerical root-finding. Second, because the first dimension has a certain degree of freedom from factor matrix rotation, $\boldsymbol{V}$ can be chosen in favor of the numerical integrations that follow (see § 4.2). The integration on the first factor is precise and computationally inexpensive, regardless of the choice of $\boldsymbol{V}$, because it is analytic. Last, analytic pricing can capture the tail probability (e.g., option price of a far-out-of-the-money strike), which Monte Carlo simulation or discretization-based numerical integration cannot easily do.

2. Quadrature integration on other dimensions

Integration over other dimensions $\dot{\boldsymbol{z}}$ can be performed using a numerical quadrature. As the integration is weighted by normal distribution density, Gauss-Hermite quadrature (GHQ) is the most suitable choice. Let $\{\dot{\boldsymbol{z}}_{m}\}$ and $\{h_{m}\}$ be the points and weights, respectively, of the GHQ associated with $n(\dot{\boldsymbol{z}})$, generated over the dimensions, $(z_{2},\cdots,z_{N})$. Subsequently, the option price becomes a weighted sum as follows:

Here, $M$ is the total number of nodes, $M=\prod_{j=2}^{N}M_{j}$, where $M_{j}$ is the node size of the $j$-th dimension. Therefore, the option price is casted on a linear combination of asset prices into the weighted sum of the multi-asset BSM prices in (9), where the forward prices of the assets vary as $F_{k}\,f_{k}(\dot{\boldsymbol{z}}_{m})$.

Moreover, the same integration scheme can be applied to the quantities of interest other than the call option price. A few examples are given below.

3. Dimensionality reduction

Since the problem involves the factor matrix $\boldsymbol{V}$, it naturally leads to the possible dimensionality reduction of low varying factors in the context of PCA. To minimize the side effects, the following approach is adopted. Without loss of generality, assume that the factor strengths $|\boldsymbol{V}_{*j}|$ are in decreasing order of $j$ and the aim is to reduce the dimensions of $N^{\prime}<j\leq N$ because these $|\boldsymbol{V}_{*j}|$ are small. Thus, it is assumed that the variation of $d(\dot{\boldsymbol{z}})$ on such $z_{j}$ is also small and the dependence on these dimensions is ignored as $d(\dot{\boldsymbol{z}})\approx d(\ddot{\boldsymbol{z}})$, where $\ddot{\boldsymbol{z}}=(0,\,z_{2},\cdots,z_{N^{\prime}},\,0,\cdots)^{\scriptscriptstyle T}$ is the state vector of the surviving dimensions, with zeros padded to the rest for convenience. In addition, notation $d(\ddot{\boldsymbol{z}})$ should be similarly interpreted as $d(z_{2},\cdots,z_{N^{\prime}})$. Because the dependence of $C_{\text{BS}}(\dot{\boldsymbol{z}})$ on the reduced dimensions occurs only through $f_{k}(\dot{\boldsymbol{z}})$, the integration of (11) over the truncated dimensions can be moved to (9) instead, which can be done analytically. With abuse of notation, the integral of $f_{k}(\dot{\boldsymbol{z}})$ over those dimensions can similarly be given as

Thus, the forward price can be preserved as $F_{k}=E(F_{k}\,f_{k}(\ddot{\boldsymbol{z}}))$, even after dimensionality reduction. The previous results—(8), (9), and (11)—remain remarkably consistent under reduced dimensions through pure notational changes—$N$ to $N^{\prime}$, $\dot{\boldsymbol{z}}$ to $\ddot{\boldsymbol{z}}$, and $f_{k}(\dot{\boldsymbol{z}})$ to $f_{k}(\ddot{\boldsymbol{z}})$. In particular, through quadrature integration (11) over $\ddot{\boldsymbol{z}}$, the total number of nodes is reduced to $M=\prod_{j=2}^{N^{\prime}}M_{j}$. This dimensionality reduction is critical for pricing Asian options, as discussed later in § 5.

4. Forward price as control variate

In case of too sparse quadrature nodes, the error from integration can be non-negligible. This error can be reduced by using the forward price $F_{k}$ as control variate. Let $\bar{f_{k}}$ be the numerically evaluated expectation of $f_{k}(\dot{\boldsymbol{z}})$, $\bar{f_{k}}=\sum_{m=1}^{M}h_{m}\,f_{k}(\dot{\boldsymbol{z}}_{m})$, which is not exactly equal to 1. For example, for a standard normal $z$, $E(e^{-\frac{1}{2}+z})$ deviates by $-6.3\times 10^{-3}$ under the GHQ evaluation with three nodes; it also deviates by $-4.6\times 10^{-4}$ with four nodes from the true value $1$. Thus, $F_{k}$ is mispriced by $F_{k}(\bar{f_{k}}-1)$ under the GHQ evaluation. This error is also present in the put-call parity of (11) and (12):

This leads to problems such as inconsistent implied volatilities for put and call options with the same strike price. Since the sensitivity (delta) of each $F_{k}$ is computed as (14), the option prices can be adjusted by using them as the coefficients of the control variate:

The put-call parity of the adjusted prices, $C^{\prime}$ and $P^{\prime}$, holds exactly.

Optimal choice of risk factor matrix

Through a simple example, it is demonstrated why a naive numerical integration suffers from slow convergence and how a proper rotation of the factor matrix can improve convergence. Consider a basket put option on two uncorrelated assets ($N=2$) with parameters $\boldsymbol{\Sigma}=\boldsymbol{I},r=0,w_{k}=1,F_{k}=e^{1/2}$ and $q_{k}=0$ for $k=1,2$ to ensure that the payoff is given as $(K-e^{x_{1}}-e^{x_{2}})^{+}$ for the independent standard normal RVs, $x_{1}$ and $x_{2}$. A put option is chosen to clearly illustrate the singularity from the vanishing exercise region. However, the same result holds for the call option also owing to put-call parity. The put option price is

where $P_{\text{BS}}(x_{2})$ is the integration of the payoff along the $x_{1}$ axis from $x_{1}=-\infty$ to $-d(x_{2})=\log(K-e^{x_{2}})$. As shown in Fig 1(a), the exercise boundary $-d(x_{2})$ diverges to $-\infty$ as $x_{2}$ approaches $\log(K)$ from the left; thus, $P_{\text{BS}}(x_{2})=0$ for $x_{2}\geq\log(K)$. In order to accurately evaluate the numerical integration over the $x_{2}$ axis, the discretization around the singularity at $x_{2}=\log(K)$ should be dense.

Alternatively, consider the $45^{\circ}$-rotated coordinate $(z_{1},z_{2})$ under which the payoff becomes $(K-e^{(z_{1}-z_{2})/\sqrt{2}}-e^{(z_{1}+z_{2})/\sqrt{2}})^{+}$. This yields the exercise boundary, $z_{1}=-d(z_{2})=-\sqrt{2}\log\big(2\cosh(z_{2}/\sqrt{2})/K\big)$, and the option price becomes

Since the boundary $-d(z_{2})$ exists at all $z_{2}$, $P_{\text{BS}}(z_{2})$ is infinitely differentiable in all $z_{2}$ and hence suitable for numerical integration along $z_{2}$ (see Fig 1(b) for $P_{\text{BS}}(x_{2})$ and $P_{\text{BS}}(z_{2})$.) As shown in Fig 1(c), the error from quadrature integration under $(z_{2},z_{1})$ decreases exponentially as the number of nodes increases, but the error under $(x_{2},x_{1})$ decreases very slowly.

2. Selection of first factor

The aforementioned intuition is generalized to the correlated and higher dimensional cases. The following two criteria are set for the selection of $\boldsymbol{V}$: (i) the exercise boundary $d(\dot{\boldsymbol{z}})$ of (8) should exist for all $\dot{\boldsymbol{z}}$ and $K$, and (ii) the variation of $d(\dot{\boldsymbol{z}})$, $|\partial d(\dot{\boldsymbol{z}})/\partial z_{j}|$, for $j\geq 2$, should be minimized. The purpose is to make $C_{\text{BS}}(\dot{\boldsymbol{z}})$ not only differentiable over $\dot{\boldsymbol{z}}$, without $-d(\dot{\boldsymbol{z}})$ diverging, but also low varying to the maximum extent possible.

As the coefficient functions $\{f_{k}(\dot{\boldsymbol{z}})\}$ can take almost any arbitrary positive values, $w_{k}{V}_{k1}>0$ is imposed for all $k$ as a sufficient condition to satisfy (i). In such a constraint, the left-hand side of (8) represents a strictly monotonic function of $z_{1}$, with the value range of $(0,\infty)$ for basket and Asian options or $(-\infty,\infty)$ for spread options. Hence, a unique root $d(\dot{\boldsymbol{z}})$ exists for any $\dot{\boldsymbol{z}}$ and non-trivial $K$.

For (ii), the following linearized approximation of the GBM is applied: $\exp(-\frac{1}{2}{V}_{kj}^{2}+{V}_{kj}z_{j})\approx 1+{V}_{kj}z_{j}$, assuming a small variance, $\|\boldsymbol{V}\|_{F}\ll 1$. After ignoring the second-order and higher terms, (8) approximates

and $\partial d(\dot{\boldsymbol{z}})/\partial z_{j}$ is obtained as the constant

Here, $\boldsymbol{g}$ is the normalized forward-adjusted weight vector, $g_{k}\propto w_{k}F_{k}$, where $|\boldsymbol{g}|=1$. The partial derivatives are minimized to zero when $\boldsymbol{Q}_{*1}$ is aligned to the direction of $\boldsymbol{C}^{\scriptscriptstyle T}\boldsymbol{g}$ because the choice maximizes the denominator and makes the numerator zero from the orthogonality between $\boldsymbol{Q}_{*1}$ and $\boldsymbol{Q}_{*j}$ for $j\geq 2$. Therefore, the optimal first factor is determined:

This is equivalent to rotating the $z_{1}$ axis to the steepest ascending direction in the left-hand side of (20), which is in agreement with the observation from § 4.1. The other axes span the slowly varying dimensions, which are subject to costly numerical integrations. The overall computation cost can be minimized in this manner.

However, $\boldsymbol{V}_{*1}$ in (22) does not always conform to the earlier constraint $w_{k}{V}_{k1}>0$. In the case of $w_{k}{V}_{k1}\leq 0$, for certain $k$, ${V}_{k1}$ is adjusted by pushing it into the conforming region by

for a small $\varepsilon>0$ and rescaling factor $\mu>0$, thereby making $\boldsymbol{Q}_{*1}=\boldsymbol{C}^{-1}\boldsymbol{V}_{*1}^{\text{(adj)}}$ a unit vector ($\mu=1$ if no adjustment). Here, $\sqrt{{\Sigma}_{kk}}$ is used as a characteristic scale of ${V}_{k1}$ because $|{V}_{k1}|\leq\sqrt{{\Sigma}_{kk}}$.

3. Remaining factors

The remaining columns $\boldsymbol{V}_{*j}$ for $j\geq 2$ are determined using singular value decomposition (SVD); this rearranges the columns orthogonally and in decreasing order of factor strength $|\boldsymbol{V}_{*j}|$. This process is executed in the following two steps. First, an orthonormal rotation matrix is found, whose first column is the same as $\boldsymbol{Q}_{*1}$. The computationally lightest choice is the Householder reflection matrix $\boldsymbol{R}$, which maps $\boldsymbol{e}_{1}=(1,0,\cdots)^{\scriptscriptstyle T}$ to $\boldsymbol{Q}_{*1}$ using the mirror image

Thus, the first column of $\boldsymbol{C}\boldsymbol{R}$ is equal to $\boldsymbol{V}_{*1}$. Second, the remaining columns of $\boldsymbol{C}\boldsymbol{R}$ are rearranged via the reduced-size SVD, $\boldsymbol{C}\,(\boldsymbol{R}_{*2}\;\cdots\;\boldsymbol{R}_{*N})=\dot{\boldsymbol{U}}\,\dot{\boldsymbol{D}}\,\dot{\boldsymbol{Q}}^{\scriptscriptstyle T}$, where $\dot{\boldsymbol{D}}$ is an $(N-1)\times(N-1)$ diagonal matrix with the (non-negative) singular values in decreasing order, and $\dot{\boldsymbol{U}}$ and $\dot{\boldsymbol{Q}}$ are the $N\times(N-1)$ and $(N-1)\times(N-1)$ matrices, respectively, satisfying $\dot{\boldsymbol{U}}^{\scriptscriptstyle T}\dot{\boldsymbol{U}}=\dot{\boldsymbol{Q}}^{\scriptscriptstyle T}\dot{\boldsymbol{Q}}=\boldsymbol{I}$. Finally, the full $\boldsymbol{V}$ is obtained by the following column-wise concatenation

and the corresponding $\boldsymbol{Q}$ as

where $\boldsymbol{0}$ is the zero vector. Because of the linearized assumption, the choice of $\boldsymbol{V}$ is independent of the strike price $K$; this ensures that if $\boldsymbol{V}$ is computed once, then it can be used for options with multiple strike prices.

Remarks on Asian options

This section discusses a few implications of the method employed here in the context of Asian options. First, since the covariance $\boldsymbol{\Sigma}$ is from the self-correlation in a Brownian motion, the Cholesky decomposition $\boldsymbol{C}$ is simply computed as

Therefore, there is no computational burden for a large $N$.

Second, all elements of $\boldsymbol{V}_{*1}$ in (22) are positive because all elements of $\boldsymbol{\Sigma}$ and $\boldsymbol{g}$ are positive for Asian options. Essentially, the selected $\boldsymbol{V}_{*1}$ is not compromised by the adjustment step of (23).

Third, the columns of factor matrix $\boldsymbol{V}$ can be interpreted as a series representation of the Brownian motion $\sigma W(t)$ on the discretized time set $\{t_{k}\}$, $\sigma\,\big(W(t_{1}),\cdots,W(t_{N})\big)^{\scriptscriptstyle T}\,{\mathrel{\mathop{\kern 0.0pt=}\limits^{d}}}\,\boldsymbol{V}\boldsymbol{z}$. If $V_{j}(t)$ is the continuum limit of ${V}_{kj}$, as $N\rightarrow\infty$, with $t=k/N$ fixed for $\sigma=1$, then the set $\{V_{j}(t)\}$ would serve as a series expansion of $W(t)$:

for independent standard normals $\{z_{j}\}$. Therefore, it is worth comparing this expansion to the well-known Karhunen–Loève expansions:

Figure 2 shows the first three factors from the two expansions. The terms are similar, but with a subtle difference in the first factor. For the normalized parameters, that is, for $\sigma=1$, $t_{k}=k/N\;(T=1)$, $w_{k}=1/N$, and $r=q=0$, ${V}_{1}(t)$ is given by the continuous version of (22):

Since the Karhunen–Loève expansion is the PCA of $W(t)$ in the functional space, $V_{1}^{\text{(KL)}}(t)$ is chosen to maximize the $L^{2}$-norm, $\int_{0}^{1}{V}_{1}^{2}(t)dt$. Therefore, $\|V_{1}^{\text{(KL)}}(t)\|_{2}=2/\pi\,(\approx 0.6366)$ is larger than $\|{V}_{1}(t)\|_{2}=\sqrt{2/5}\,(\approx 0.6325)$. On the other hand, $V_{1}(t)$ is chosen to maximize $\int_{0}^{1}V_{1}(t)dt$. Thus, $\int_{0}^{1}V_{1}(t)\,dt=1/\sqrt{3}\,(\approx 0.5774)$ is larger than $\int_{0}^{1}V^{(\text{KL})}_{1}(t)\,dt=4\sqrt{2}/\pi^{2}\,(\approx 0.5732)$. The continuous representation (29) is valid only under constant weight, and no longer holds if, for example, $r\neq 0$ or $w_{k}$ is not constant. Therefore, generally $\boldsymbol{V}_{*1}$ is numerically computed instead of using (29).

Fourth, dimensionality reduction is critically effective in Asian options because the dimension $N$ is large. For a one-year maturity, the monthly averaging, weekly averaging, and daily averaging correspond approximately to $N=12$, $N=50$, and $N=250$, respectively. Even if a few nodes per dimension are used, the total number of nodes becomes prohibitively large. However, the first several factors explain most of the variance in option prices, while the remaining factors have negligible impact. Table 1 shows the portion of the explained variance, $\sum_{j=1}^{N^{\prime}}|\boldsymbol{V}_{*j}|^{2}/\|\boldsymbol{V}\|_{F}^{2}$, as a function of increasing $N^{\prime}$. The first factor, $\boldsymbol{V}_{*1}$, accounts for the largest part, about 80%, with the cumulative portion reaching about 96% at $N^{\prime}=5$. The numerical results in § 6 show that five dimensions, one under analytic and four under numerical integration, are sufficient for accurate option pricing. If the averaging feature is meant to avoid market manipulation, as is often the case, averaging should start sometime near maturity. Since the correlation in this case is overall higher compared to that in the case of immediate averaging, the explained variance portion is higher and the dimensionality reduction more effective.

Finally, the valuation of continuously monitored Asian options is cast into discrete monitoring. For a better integration convergence over time, the Simpson’s rule weights are used rather than constant weights. For the discretization step $\Delta T$ such that $N=T/\Delta T$ is even, the observation time and weights are given as

with the exception of the weights at the two end points, $w_{0}=w_{N}=\Delta T/(3T)$.

Numerical results

The method is implemented in R (Ver. 3.3.2, 64-bit) on a personal computer running the Windows 10 operating system with an Intel core i5 2.2 GHz CPU and 8 GB RAM. Seven parameter sets are tested, the first six of which are described in Tables 2 and 3. The sets are labeled as S for spread, B for basket, and A for Asian options. The set A3 for continuously monitored Asian options is separately displayed in Table 12 along with the results. Except for set S2, the parameters are the same ones used in previous studies to enable easier comparison.

The numerical results are reported in Tables 4$\sim$12. Except for the three Asian option sets, two versions of the prices are computed. One version is the “fast” price, for which the minimal GHQ nodes are used with a target precision of 3$\sim$4 decimals, which is both practical and sufficient. The computational cost for the fast price is inexpensive. The other version is the converged price; the fast price error is measured against this. The converged price is obtained as the node sizes are increased. Seven decimal precisions are targeted for the converged price. For the Asian option cases (Tables 10$\sim$12), only the fast price is reported; the error is measured from the previous studies: Černỳ and Kyriakou 2011 for the discrete cases (A1, and A2) and Linetsky 2004 for the continuous case (A3). All prices in the tables are the present values of the call options; that is, the forward value discounted by $e^{-rT}$.

Some tables show the used factor matrix $\boldsymbol{V}$ in the following representation to provide extra properties besides the matrix itself.

The center itself is $\boldsymbol{V}$. The upper and right-hand side panels show the $L^{2}$-norm of the columns and rows, respectively, while the upper right corner is the Frobenius norm. The left-hand side panel shows the forward-adjusted weight vector $\boldsymbol{g}$, and the upper left corner shows the dot product $\boldsymbol{g}^{\scriptscriptstyle T}\boldsymbol{V}_{*1}$. The lower panel shows the node size $M_{j}$ for the $j$-th dimension of the $j\geq 2$, with the total size $M$ in the lower right corner.

For implementation, $M_{j}$ must be chosen in an economic manner. The node sizes need not be the same for all dimensions. Therefore, $M_{j}$ can be selected to ensure that all dimensions homogeneously achieve a similar level of accuracy. A general guideline is that the density of the nodes should be proportional to $|\partial d(\dot{\boldsymbol{z}})/\partial z_{j}|$ in order to effectively capture the change in $d(\dot{\boldsymbol{z}})$. To this extent, the following rule is used to systematically determine $M_{j}$ in the numerical tests:

where $[x]$ is the nearest integer of $x$ and $\lambda$ is the coefficient for the level of accuracy. The ratio $|\boldsymbol{V}_{*j}|/|\boldsymbol{g}^{\scriptscriptstyle T}\boldsymbol{V}_{*1}|$ is obtained by applying the Cauchy-Schwarz inequality, $|\boldsymbol{g}^{\scriptscriptstyle T}\boldsymbol{V}_{*j}|<|\boldsymbol{g}|\cdot|\boldsymbol{V}_{*j}|=|\boldsymbol{V}_{*j}|$, to (21), which is thereby understood as an approximate upper bound of $|\partial d(\dot{\boldsymbol{z}})/\partial z_{j}|$. This rule serves also as a criterion for dimension reduction: if $M_{j}=1$ for some $j$, the dimension can be truncated according to § 3.3. Moreover, this rule is independent of $K$; it ensures that if $\{\dot{\boldsymbol{z}}_{m}\}$ and $\{w_{m}\}$ are computed once, they can be used for options with multiple values of $K$. For Asian options, however, $M_{j}=3$ for $2\leq j\leq 5$ ($M=81$) is empirically found to work very well; thus, the test cases for Asian options do not resort to (32).

In S1, the spread call options are priced for varying strikes from $K=0$ (in-the-money) to $4$ (at-the-money). For the speed of convergence, the node size is increased for the second dimension from $M_{2}=2$. As shown in Table 4, convergence is extremely fast. While the prices with $M_{2}=2$ are already accurate, they converge within seven decimals at $M_{2}=3$. The table also shows that the control variate correction of § 3.4 can further reduce error. While Hurd and Zhou 2010 and Caldana and Fusai 2013 show accuracy similar to the $M_{2}=3$ result, their methods have to evaluate the expensive Fourier inversion in two and one dimensions, respectively. The risk factor matrix is presented in Table 4(b). In this example, the adjustment step (23) is triggered for ${V}_{21}$ with $\varepsilon=0.01$. In addition, this study independently implements the analytic approximation methods of Bjerksund and Stensland 2014, Lo 2015, and Li et al. 2008 for spread options. For Lo 2015, the “Strang’s splitting approximation I” method in the reference is implemented; this is given in an explicit formula. All the three methods work well on S1, showing errors in the order of $10^{-6}$ or less.

In S2, the at-the-money spread call options are priced for varying correlations $\rho_{12}$ from $90\%$ to $-90\%$. The results are shown in Table 5(a). This parameter set is designed to make the first asset to follow a displaced GBM, $dS_{1}(t)=\sigma_{1}(S_{1}(t)+L)\,dW_{1}(t)$ with $L=100$ to demonstrate a negative implied volatility skew. Therefore, the true parameters are similar to those of $S_{2}$, as $S_{1}(0)=100$ and $\sigma_{1}\approx 30\%$, and $K=100$ corresponds to the at-the-money strike price. The factor matrix $\boldsymbol{V}$ changes with the change in $\rho_{12}$ and $M_{2}$ is determined according to (32). For a fast price, $M_{2}$ varies from $17$ to $2$ when $\lambda=3$ is used. Since the components in $\boldsymbol{V}_{*1}$ must have opposite signs for spread options, its norm weakens under positive correlation, thus requiring denser nodes. The quadrature size rule (32) works reasonably well, although it tends to over-allocate nodes under high correlation. The price quickly converges as $\lambda$ is increased, to achieve the seven-decimal precision at $\lambda=9$. Tables 5(b) and (c) show $\boldsymbol{V}$ for $\rho_{12}=90\%$ and $-90\%$, respectively. Note that the two columns switch places with a sign change. The performance of analytic approximation methods in S2 is not as good as that in S1; this might be attributed to the displaced GBM features, that is, the significant difference between two spot prices and the big strike price. Although the method of Li et al. 2008 is much more accurate than the other two methods, it illustrates the limitation of analytic approximation.

2. Basket Option

Parameter sets B1 and B2 have been frequently used as benchmarks in previous studies. Here, the convergent option values are reported for the first time. Parameter set B1 is extracted from Krekel et al. 2004, where the performance results of several analytic approximation methods are compared. Tables 6 and 7 provide the results for varying $K$ and $\rho_{k\neq j}$, respectively, from the base parameter values. Moreover, Table 8 tests the inhomogeneous volatilities by varying $\sigma_{k}$ simultaneously for $1\leq k\leq 3$, while keeping $\sigma_{4}=100\%$ fixed. In all tests, $\lambda=9$ give consistent fast prices, which are more accurate than the benchmark Monte Carlo prices of Krekel et al. 2004; moreover, the computation cost is much cheaper. The test in Table 6 shows that $\lambda=9$ corresponds to $M_{j}=5$ for $2\leq j\leq 4$ ($M=125$). In addition, the node sizes in Tables 7 and 8 are shown in their respective last columns. Contrary to the spread option, a lower correlation requires denser nodes in a basket option. In the test, the lower bound of the correlation $\rho_{k\neq j}$ is $-1/3$; this is to ensure that the covariance matrix is positive-semidefinite. Also note that only a few nodes per dimension produce accurate prices in the inhomogeneous volatility test. No approximation method produces satisfactory precision in the same test (Krekel et al. 2004).

The parameter set B2, referred to as the G-7 indices basket option, has been tested in Milevsky and Posner 1998b and Zhou and Wang 2008. The call option prices for varying $K$ and $T$ are shown in Table 9(a) and the factor matrix $\boldsymbol{V}$ is shown in Table 9(b). The current method works well for this $N=7$ case. The error of the fast prices with $\lambda=3$ ($M=432$) is in the order of $10^{-4}$ at most, which is smaller than the standard error of the Monte Carlo simulation in Zhou and Wang 2008. The seven-digit convergent prices are computed with $\lambda=12$ ($M\approx 1.2\times 10^{5}$).

3. Asian options

Three parameter sets for the Asian options are tested. The results for the discrete monitoring sets A1 and A2 are reported in Tables 10 and 11, respectively. A parameter set, which is almost the same as A1, except for $r=4\%$, has been popularly tested in the literature as well; for example, Černỳ and Kyriakou 2011; Fusai et al. 2011 and Cai et al. 2013. They should not be confused. In Table 10, $K$ and $\sigma$ are varied, whereas in Table 11, $K$ and $N$ are varied. The errors are measured from Černỳ and Kyriakou 2011, where the results are reported with an accuracy of $10^{-7}$. The fast prices computed with only 81 nodes over five dimensions ($M_{j}=3$ for $2\leq j\leq 5$) have errors in the order of $10^{-4}$ or less. The overall underpricing is due to a decline in variance from dimensionality reduction. In A1 and A2, the average computation time per option price is $0.01$, $0.02$, and $0.06$ seconds for $N=12,\;50$, and $250$, respectively. This can be accelerated when the prices for multiple values of $K$ are computed together because the computations for $\boldsymbol{V}$, $\{\dot{\boldsymbol{z}}_{m}\}$, and $\{h_{m}\}$ are not repeated. Although a direct comparison on CPU time is difficult to obtain because of difference in computing environments, Černỳ and Kyriakou 2011 report 1 to 0.3 seconds as $\sigma$ varies from $10\%$ to $50\%$ for computing $N=50$ cases with a five-decimal precision. Table 12 reports the result for the continuous monitoring set A3. Because time is discretized with $\Delta T=1/200$, $N=200$ for $T=1$ and $N=400$ for $T=2$. The fast prices computed with 81 nodes match the values of Linetsky 2004, computed with up to 10-decimal accuracy. The computation time per option is $0.04$ and $0.18$ seconds for $N=200$ and $400$ cases, respectively.

Conclusion

Option pricing under a multivariate BSM model is a challenging task because of the curse of dimensionality, and has a long history of research. This study eases this curse significantly by replicating the option prices as a weighted sum of single-factor BSM prices under an optimally rotated state space. Moreover, this method can be uniformly applied to spread, basket, and Asian options. Numerical examples show that this method is both fast and accurate.

References