Primal-dual active-set method for solving the unilateral pricing problem of American better-of options on two assets

1.   Introduction

Due to the increased awareness of risk aversion, the financial derivatives (e.g., options) have attracted more and more attention. An option is a derivative, which can be exercised on the maturity date (European option) or at any time prior to the maturity date (American option). Based on the geometric Brownian motion, the well-known Black-Scholes equation (BS equation) for pricing options was established in 1973. Whereafter, a variety of closed-form expressions have been derived for most European options ^[1,2,3]. Up to now, there has been no similar results for American options yet. By virtue of the flexibility in the choice of the exercise time, American options have been developed in the option market rapidly. Originally, traders and researchers focused on single-asset options, which are of simple structures ^[4]. With the development of the financial market and the increasing demand of investors, plentiful multi-asset options have emerged, and the corresponding American options were studied intensively. Although multi-asset American options can be divided into some categories, and may cover any number of assets, the better-of option on two assets, which belongs to the rainbow option family, is typical and popular because of its simplicity and practicality.

The better-of option was first named "option on the maximum of two risky assets" by Stulz in 1982 ^[5], and then was called its present name by Jiang ^[6]. By exploiting the relationship of the pricing formulas for two single asset options, Stulz presented the closed-form solution for the European better-of option. Jiang stated that the American better-of option satisfies a multi-dimensional BS equation with a special payment function. In most instances, the underlying assets of the options pay dividends or have other cash outflows. As is well known, the standard American option written on a single asset that pays dividends can be exercised optimally before the maturity date ^[7]. And so are multi-asset dividend paying options ^[8]. When the underlying assets pay continuous dividend yields, the set of optimal prices of multi-asset options with respect to the time forms a continuous surface, which is called the optimal exercise boundary. In this paper, we consider a special two-asset American better-of option pricing problem, where only one underlying asset pays dividends. The optimal exercise boundary for the concerned model divides the solving domain into two parts: one is the bounded region called the holding domain, and the other is the unbounded one called the exercising domain.

In practice, various methods have been used to price options accurately and efficiently. They can be classified into two categories: analytical approximations ^[9,10] and numerical solutions ^[11]. On account of no analytical formula, the pricing problem for American options has to be solved numerically. Binomial method (BM), Monte Carlo (MC) simulations, finite difference methods (FDM) and finite element methods (FEM) ^[12] are commonly used numerical methods for pricing American options. BM is the most classic numerical technique for options. In 1979, Cox, Ross and Rubinstein first applied the binomial model to price American options ^[13]. Five years later, Amin and Khanna proved the convergence of this method ^[14]. With the increase in the number of the underlying assets, the high computational cost and the slow convergence rate incurred by BM. Based on the pioneering works of Boyle, Bossaerts and Tilley, MC simulations are frequently used to the American option pricing problems ^[15]. The advantage of MC is that it does not depend on the dimension of the problem and thus does not suffer from the curse of dimensionality. As a generalization of BM, FDM is introduced to solve the pricing problems and is widely used because of its simple form. The FDM is easy for implementation and efficient for problems with regular domain. While for problems with irregular domain, FEM is more flexible. It can handle the irregular domain as well as complex boundary conditions efficiently. Meanwhile, FEM has a solid theoretical framework with robust numerical performance. Therefore, FEM has became more and more popular in classical option pricing problem ^[16], and has been successfully used to deal with better-of options recently ^[17,18]. In this paper, we shall adopt FEM to discretize the two-asset American better-of option pricing problem when only one asset pays dividend.

The concerned two-asset American better-of option in this paper satisfies a two-dimension parabolic linear complementary problem (LCP) on an unbounded domain. To solve this model efficiently, we must think about the following issues: 1) In order to better characterize the singularity of the option price near the maturity date, the design of the spatial mesh must be considered carefully. 2) Due to the solution domain being unbounded, it is difficult to design algorithms directly. The truncation method and the corresponding boundary condition, which can guarantee the accuracy of the solution, are important. 3) For the dependency relationship of the option price and the exercise boundary, the pricing model becomes a highly nonlinear problem. The efficiency of numerical algorithm used to settle the option price and the optimal exercise boundary simultaneously is the key issue.

As regards the first difficulty, the geometric grid for dealing with the classical single asset American options is adopted to guarantee the accuracy of the better-of option around the singularity. For the second issue, there are two mainstream techniques: the transparent boundary condition method and the far-field boundary condition method. The former one was proposed to solve classical American options by Han and Wu in 2003 ^[19], and improved by Ehrhardt and Mickens in 2008 ^[20], which is suitable for use in combination with FDM. The latter method was proposed by Kangro and Nicolaides ^[21]. They presented a reasonable boundary condition along the truncation boundary, and established the pointwise estimate of this truncation method via the comparison principle. Whereafter, many works for pricing American options using this method have been published ^[16,18]. We shall follow this approach to deal with our problem. Regarding the third issue, the primal-dual active-set (PDAS) method shall be used to solve it efficiently. The PDAS method is a special case of the generalized Moreau-Yosida approximations ^[22], and is also proved to be a special case of the Newton-type method under some moderate conditions ^[23]. It is efficient for quadratic programming problems and LCPs. Moreover, the global convergence of PDAS method has been proved by Bergounioux et al. in ^[24]. For the applications of PDAS on option pricing problems, we refer to ^[17,25] and references therein for the rich literature.

The concerned model in this paper is the same as the problem in ^[18], but the numerical algorithm is very different. Although both methods can be used to obtain the option price and the optimal exercise boundary simultaneously. Compared with the Newton's method in ^[18], our proposed algorithm has obvious advantages in both theoretical analysis and numerical performance. Firstly, the convergence of our method can be analysed, while the Newton's method in ^[18] can not be. Secondly, PDAS method as an active-set strategy could effectively reduce the number of the iterations in each time layer, and leads to rapid convergence, which can be verified in the numerical simulations. Moreover, PDAS method is faster than the projected successive overrelaxation method. Therefore, the proposed algorithm in this paper is an efficient method for pricing two-asset American better-of options when only one asset pays dividends.

The arrangement of this paper is as follows. In Section 2, the variational inequality (VI) model for American better-of options is introduced firstly. Based on some conventional numeraire transformations and far-field technique, the bounded LCP corresponding to the VI is established. In Section 3, the full discretization model of the LCP is constructed by FDM with uniform partition in temporal direction and FEM with geometric mesh in spatial direction, respectively. Then, the PDAS method is adopted to solve the resulting discretized optimization system. The convergence analysis and numerical experiments are implemented in Section 4. The concluding remarks are given in Section 5.

2.   The pricing models

In this section, we shall introduce the VI model for American better-of options on two underlying assets. Furthermore, based on some changes of variables and the far-field truncation technique, the associated LCP shall be presented.

2.1.   The variational inequality model

For the simplicity of the expression, we introduce some notations firstly. $S_{i}$, $\sigma_{i}$ and $q_{i} \; (i = 1, 2)$ stand for the price, the volatility, and the dividend of the $i$-th underlying asset, respectively. The correlation coefficient of these two assets is denoted by $\rho$. It needs to point out that we only consider the case where only one underlying asset pays dividend, as the case in which both underlying assets pay dividends have been discussed exhaustively ^[17]. Let $r$, $t$, and $T$ be the interest rate, the arbitrary time, and the maturity date, respectively. Assume that there is no transaction fee to pay and no arbitrage on the market. Then, the VI model corresponding to the American better-of option price $P = P(S_{1}, S_{2}, t)$ with $q_{1} = 0$ and $q_{2} > 0$ can be described as

where the payoff function $f(S_{1}, S_{2}) = \max(S_{1}, S_{2})$, the solution domain

and the operator

here, $\boldsymbol{ R}^{2}_{+} = \{(S_{1}, S_{2})\mid S_{i}\in(0, +\infty), i = 1, 2\}.$

The solution domain $\Sigma$ of the pricing model (2.1) can be divided into two parts by the optimal exercise boundary, where the payoff of an option holder reaches the maximum. Since $q_1 = 0$ and $q_2 > 0$, the exercising domain denoted by $\Sigma_2$ is on the left of the optimal exercise boundary. When the prices of the underlying assets belong to $\Sigma_2$, the price of the option is equal to the payoff value by exercising the option, i.e., $P(S_1, S_2, t) = f(S_1, S_2)$. Otherwise, the holders may suffer a loss. On the other hand, the holding domain denoted by $\Sigma_1$ is on the right of the optimal exercise boundary. On this subdomain, the price of the option is unknown, and satisfies $P(S_1, S_2, t) > f(S_1, S_2)$. Therefore, the model (2.1) can be rewritten as

It is worth to notice that there is an unknown interface between $\Sigma_{1}$ and $\Sigma_{2}$. In this paper, we denote it by $\Gamma$. Hence, the pricing model (2.1) is equivalent to a two-dimensional free boundary problem.

By means of traditional numeraire transformations

the pricing model (2.1) can be transformed into the following one-dimensional VI model

where the solution domain after the transformations

and the simplified operator

Similarly, $\hat{\Sigma}_{1}$, $\hat{\Sigma}_{2}$, and $\hat{\Gamma}$ represent the corresponding holding domain, the exercising domain, and the optimal exercise boundary, respectively. The one-dimensional VI model (2.4) can be reformulated as

Here, $\hat{\Sigma}_{1} = \big\{(s, t)\mid s\in\big(\gamma(t), +\infty\big), \; t\in[0, T)\big\}$, $\hat{\Sigma}_{2} = \hat{\Sigma}\backslash\hat{\Sigma} _{1}$ and $\hat{\Gamma} = \gamma(t)$ is the corresponding unknown free boundary. Based on the above discussions, we shall only solve the following unilateral free boundary problem

Furthermore, let

then the free boundary problem model (2.6) shall be transformed into a standard American put option pricing problem

where the optimal exercise boundary $\gamma(t)$ is a bounded and nondecreasing function ^[26], and satisfies

here, $\beta_{i} = \dfrac{\alpha_{i}-1}{\alpha_{i}}\; (i = 1, 2)$, $\alpha_{1} < 0$, and $\alpha_{2} > 0$ are the solutions of the following equation

2.2.   The bounded linear complementary problem

In order to solve the free boundary problem model (2.8) numerically, the unbounded domain must be truncated firstly. Following the ideas of Holmes and Yang ^[16], we shall introduce a far-field boundary method to deal with this problem, by which accurate solutions can be obtained efficiently.

Lemma 1. For a given positive number $\varepsilon \in (0, 1)$, let $\psi = \dfrac{1}{2}\hat{\sigma}^{2}$, $\alpha_{0} = \dfrac{r-q}{\hat{\sigma}^{2}}-\dfrac{1}{2}$, then we have

where

Proof. By applying the variable substitutions

the free boundary problem model (2.8) can be transformed into

where the coefficients $\psi = \dfrac{1}{2}\hat{\sigma}^{2}$, $\nu = 2\alpha\psi-q_{2}$, and $\mu = (\alpha + 1)q_{2}-\psi\alpha^{2} -\beta$, and the functions

From the property of $\gamma(t)$, we can infer that $b(\eta)$ is a bounded and nonincreasing function with

Taking $\alpha = \alpha_{0}$ and $\beta = \psi\alpha^{2}_{0}+q_{2}$, yields $\nu = \mu = 0$. Therefore, the pricing problem model (2.8) becomes a heat equation. By virtue of the property of $b(\eta)$ in model (2.12), the first equality of model (2.8) holds when $z > 0$. Consequently, by the Theorem 19.3.2 in ^[27], we obtain

The pricing problem model (2.8) is equivalent to a standard American put option pricing problem with the strike price $K = 1$. Based on the fact that the value of the American put option is always less than or equal to the strike price and the transformations model (2.10), we can get $w(0, \eta)\leq e^{\beta\eta}$, which yields

Furthermore, by the inverse transformations of model (2.10), we have

Finally, from the definition of $\hat{Z} = -2\psi T\alpha_{0} +2\sqrt{\alpha^{2}_{0}\psi^{2}T^{2}-\psi T ln\varepsilon}$, it is easy to verify

which implies the conclusion of this Lemma.

Now, using the error bound of Lemma 1, we truncate the solution domain. For a fixed $\varepsilon \in (0, 1)$, let

then the free boundary problem model (2.8) can be approximated by the following LCP

with constraints

It needs to point out that the choice of $L$ should ensure $e^{-L}\leq\gamma(t)$ and $e^L\geq e^{\hat{Z}}$, which makes sure that the left boundary condition of model (2.14) is accurate, and the right boundary condition is reasonable.

The LCP model (2.14) is a backward variable coefficients problem on a bounded domain. To solve it efficiently, we change it to a forward constant coefficient problem by the transformations

with $a = \dfrac{2 q_{2}-\hat{\sigma}^{2}}{2 \hat{\sigma}^{2}}$, $b = \dfrac{1}{2}\hat{\sigma}^2 a^{2}+q_2$. The specific form of the bounded LCP is as follows

with constraints

Here, $q(x, \tau) = e^{ax+b\tau}f(1-e^{x}, 0)$ and $\widetilde{T} = \dfrac{\hat{\sigma}^2}{2}T$. It is easy to get that the relationship of the optimal exercise boundaries between the linear complementary problems (2.14) and (2.16) is

So far, we have transformed the original variable coefficients pricing model on a two-dimensional unbounded domain into a BLCP on a one-dimensional bounded domain, on which we can design numerical algorithms directly.

3.   The numerical algorithm

In this section, the variational problem associated with the BLCP model (2.16) is presented firstly. Then, the resulted problem shall be discretized by the FDM and the FEM in temporal direction and spatial direction, respectively. Based on the properties of the discrete system, the PDAS method is proposed to obtain the option price and the optimal exercise boundary simultaneously. The error estimates of the proposed approach are also given at the end of the section.

3.1.   The variational problem

Before presenting the variational problem corresponding to the BLCP model (2.16), we ought to introduce some function spaces and function sets for subsequent applications. $L^{2}(I)$ stands for the space of square integrable functions on $I = [-L, L]$, and define

Lemma 2. (cf. ^[28]) If $v \in L^{2}\big([0, \widetilde{T}]; H^{2}(I)\big)$and $v_{\tau}\in L^{2}\big([0, \widetilde{T}]; L^{2}(I)\big)$, thenv is the solution of the $\rm BLCP$ model (2.16) if and only ifv is the solution of the following variational inequality

$\rm(VI)$ Find $v(\cdot, \tau)\in H^{1}_{\tau}(I)$, such that$v(x, 0) = q(x, 0)$ and

In order to apply FDM and FEM to the VI model (3.1), some notations must be introduced in advance. The spatial partition of $I = [-L, L]$ and the temporal partition of $J = [0, \widetilde{T}]$ are defined as follows:

Let $I_i: = (x_{i-1}, x_i)$ denote the spatial element and $h_i: = x_i-x_{i-1}, i = 1, \dots, N$ represent the interval length. The grid size of the spatial partition is denoted by $h: = \max\limits_{1\leq i\leq N}\; h_{ i}$. In a similar way, for each temporal element $J_j: = (\tau_{j-1}, \tau_j), \triangle\tau_j: = \tau_j-\tau_{j-1}, j = 1, \dots, M$, and $\triangle\tau: = \max\limits_{1\leq j\leq M}\; \triangle\tau_{ j}$ represent the local and overall step size, respectively.

What we should pay attention to is that options are traded more frequently near the optimal exercise boundary $\Gamma$. As regards this issue, the most desirable points for the option price are around the optimal exercise boundary with $S_{1} = S_{2}$. In other words, we ought to set most nodes near the point $s = 1$ or $x = 0$ by the numeraire transformations models (2.3) and (2.15). Therefore, in the spatial direction, we shall use a geometric grid partition and an even number of intervals $N$, which ensure the nodes should be symmetrical about $x_{\frac{N}{2}} = 0$. And we resort an isometric subdivision in temporal direction. That is to say,

3.2.   The finite element method

The main goal of this subsection is to present the discretization scheme of the VI model (3.1). First, we define the set of piecewise linear functions as follows

where $\mathcal{P}_{1}$ represents the set of polynomials of degree less than or equal to $1$. Thus, the semi-discretized approximation of the VI model (3.1) by FEM can be described as

$\rm (SDA)$ Find $v_{h}(x, \tau)\in S^{1}_{\tau}(I)$, such that $v_{h}(x, 0) = q_{I}(x, 0)$ and

Here, $q_{I}(x, 0)$ stands for the piecewise linear interpolation of $q(x, 0)$ in $S^{1}_{0}(I)$. Next, the backward Euler method is applied to the SDA model (3.3) at a fixed $\tau = \tau_{j}, j = 1, \dots, M$, the full-discretized approximation of VI model (3.1) is presented as

$\rm (FDA)$ Find $v^{j}_{h}(x)\in S^{1}_{\tau_{j}}(I)$, such that $v^{0}_{h}(x) = q_{I}(x, 0)$ and

where $v^{j}_{h}(x) = v_{h}(x, \tau_j), j = 1, \dots, M$, the operator $\partial$ and the bilinear function $b(\cdot, \cdot)$ are defined as follows

For the convenience of algorithm design, we shall derive the matrix-vector form of the FDA model (3.4). Suppose that the basis function set of $S^{1}_{\tau}(I)$ is $\{\phi_{0}, \phi_{1}, \dots, \phi_{N}\}$, and the specific representations of the basis are

Then, the finite element solutions can be reformulated as

By substituting the variables $v^j_h$ in the FDA model (3.4) by the above formula, we obtain the following inequality

where $j = 1, \dots, M$,

and

Let ${{\mathit{\boldsymbol{C}}}} = \triangle\tau{\bm \varPsi}+{\bm \varPhi}$ and ${{\mathit{\boldsymbol{D}}}}^j = {\bm \varPhi}{{\mathit{\boldsymbol{V}}}}^{j-1}-{{\mathit{\boldsymbol{H}}}}^j$, therefore the inequality model (3.5) can be simplified as

In the above matrix-vector inequality model (3.6), we can verify that the matrix ${{\mathit{\boldsymbol{C}}}}$ is positive definite when $\dfrac{h^{2}}{\triangle\tau}$ is small enough.

3.3.   The primal-dual active-set method

In this subsection, to get the option price and the optimal exercise boundary simultaneously, we adopt PDAS method to solve the MVIP model (3.6). The detailed algorithm is presented in Algorithm 1:

For $j = 1, \dots, M$, when the numerical solution ${{\mathit{\boldsymbol{V}}}}^j$ is obtained, we can get $v_h^j(x)$. By means of the transformations models (2.15), (2.7), and ($\ref{eq:reduc trans}$), we can calculate the approximation of $W(s, t_j)$, $p(s, t_j)$, and $P(S_1, S_2, t_j)$. Likewise, we obtain the optimal exercise boundary $\gamma_h{(t)}$ at $t = t_j$. We present the whole process in Algorithm 2 as:

3.4.   The error estimation

The error of our proposed method mainly comes twofold: $(a)$ The error between the original problem model (2.1) and the BLCP model (2.16), which has been presented in Lemma 1. $(b)$ The error of the BLCP model (2.16) and the MVIP model (3.6). In this subsection, we focus on the error estimate of $(b)$. In order to obtain the result of estimation, we introduce some lemmas firstly.

Lemma 3. (cf. ^[29]) If $u \in H^{2}(I)$ and $I_{h}u$ is thepiecewise linear interpolation of $u$ on the grid $I_{x}$, then wehave

For any $j = 0, 1, \dots, M$, let $e_1^{j} = v^{j}-v^{j}_h$, then the following results hold

Lemma 4. (cf. ^[30]) Suppose that$v\in L^{\infty}\big(J, H^{2}(I)\big)$, $\dfrac{\partial v}{\partial\tau}\in L^{2}\big(J, H^{1}(I)\big)$, and $\dfrac{\partial v}{\partial\tau}\in L^{2}\big(J, L^{\infty}(I)\big)$.For any $j = 1, \dots, M$, let $e_2^{j} = v^{j}-I_{h}v^{j}$, then we have

Let $v_{h, \tau}$ represent the linear interpolation of $v^j_{h}$ in the temporal direction, based on the Lemmas 3 and 4, we establish the following result.

Theorem 1. Under the assumptions of Lemma 4, if in addition, $\dfrac{\partial v}{\partial\tau}$, $\dfrac{\partial^2v}{\partial{\tau}^2} \; \in L^{\infty}\big(J, H^{1}(I)\big)$, thenwe have the error estimate as follows

Proof. For $j = 1, \dots, M$, we can calculate the following equality firstly

Let $u = v_h^j$ at $\tau = \tau_j$ in VI model (3.1) and $u = I_hv^j$ in FDA model (3.4), then we obtain

Next, we add the two inequalities in model (3.8) to both sides of model (3.7), we have

For $j = 1, \dots, M$, define

then, we can simplify the inequality model (3.9) as below

Thus, we can get

where $E_l, l = 1, \dots, 4$ are defined in Lemma 4. Therefore, the first part of the above inequality can be estimated by

Based on the inequalities (3.10) and (3.11), we have

By virtue of Poincare inequality and $e_1^j\in H^1_0{(I)}$, the following result holds

Moreover, by exploiting the complex trapezoidal formula, we can conclude that

which leads to the conclusion directly ^[31].

4.   Numerical experiments

In this section, numerical experiments are performed to illustrate the advantages of our proposed algorithm. In order to show the superiority of our algorithm in computational accuracy and speed, we shall compare it with the BM and the Newton's method (NM) ^[18].

Let us consider a one-year ($T = 1$) American better-of option with $r = 0.02$. Without loss of generality, the parameters of the first underlying asset are set to be $q_1 = 0$, $\sigma_1 = 0.3/0.4$, the parameters of the other asset are set to be $q_2 = 0.02/0.03$, $\sigma_2 = 0.4$, the correlation coefficient of these two assets is set to be $\rho = 0.6/0.7$. Suppose the truncation accuracy $\varepsilon = 10^{-6}$ in model (2.14) via the definition model (2.13), the estimate model (2.9), and Lemma 1, we obtain a table about the truncation length in deferent cases.

For convenience, we choose $L = 2.8$ in all experiments, which can ensure the truncated domain cover all the solution domains corresponding to Table 1. Because there are no closed-form solution for American better-of options, we shall apply the solutions obtained by BM with 100,000 points in temporal direction as the benchmarks. The parameters in PDAS method are set to be ${\bm \lambda} = {{\mathit{\boldsymbol{0}}}}$, $\delta = 10^6$, $\varepsilon_V = 10^{-6}$. All the experiments are implemented by Matlab R2014a and on a computer with Intel Core i5 CPU of 2.2GHz.

4.1.   The optimal exercise boundary

In this subsection, we mainly verify the superiority of PDAS method in computing the optimal exercise boundaries. The error is measured by the ${L}^{2}$-norm of $\gamma_h(t)-\gamma(t)$, where $\gamma_h(t)$ and $\gamma(t)$ are the numerical solutions and the benchmark solutions, respectively. With $\rho = 0.6$ and $0.7$, Tables 2 and 3 compare the accuracies achieved by running BM, NM, and PDAS method for about the same amount of time ($\pm1\% -\pm 0.1\%$). Similarly, Tables 4 and 5 list the computation costs of the using above three methods in different cases.

From the results listed in Tables 2–5, we can conclude that $(a)$ with the similar computing time ($\pm 1\%$), the errors computed by PDAS method are one order of magnitude smaller than BM and NM. $(b)$ with the similar computing time around ($\pm 0.1\%$), the computational efficiency of PDAS is almost $100$ and $10$ times of BM and NM, respectively. Figure 1 shows the optimal exercise boundaries solved by PDAS and the benchmark BM. It is easy to find out that the numerical solutions approximate the benchmarks very well. All of the results confirm that PDAS method is an efficient algorithm to obtain the optimal exercise boundary. At the end of this subsection, to give the investors some intuitive and valuable guidance, we present the original optimal exercise surfaces $\Gamma$ of model (2.2) in Figure 2 with different parameters, respectively.

4.2.   The option price

In this subsection, we illustrate the efficiency of the PDAS method for solving option prices. Similar to the analysis in previous subsection, we compare the error estimates and computation costs in Tables 6–9.

The results about option prices listed in Tables 6–9 can fully illustrate the advantages of PDAS method in computational accuracy and speed. In order to get a more intuitive sense of option pricing, some images of the option price $p(s, 0)$ after the transformation model (2.3) and the original option price $P(S_1, S_2, 0)$ are shown in Figures 3 and 4, which are the main concerns for financial institutions. From the results in Figure 3, we confirm that the solutions of our proposed method approximate the benchmarks very well, and are all larger than the payoff $\max \{s, 1\}$, which further illustrate the efficiency and rationality of our proposed PDAS method.

4.3.   The convergence order

The practicability of PDAS method has been checked, now we shall verify the error estimates obtained in Theorem 1. Taking the temporal direction as an example, suppose that the convergence order in the temporal direction is $r$, then the logarithm of the error between the true and numerical solutions is linear with respect to the logarithm of $M$, which is the degree of freedom in the temporal direction. Therefore, we can conclude that the slope of the logarithm of the error with respect to $\log M$ is equal to $-r$. According to the result in Theorem 1, we need to verify $r = 1$ in the temporal and spatial direction, respectively. So as to reveal the relationship between the numbers $N$ of the spatial nodes and the spatial error, we choose $M = 200$ in temporal direction, which can ensure that the spatial error is dominant. Likewise, to verify the relationship between the degrees of freedom $M$ in temporal direction and the temporal error, we choose $N = 2500$. We also show blue line in each figure as the standard lines to verify the conclusion of theorem. The convergence results of options in spatial and temporal directions are shown in Figures 5 and 6.

It needs to be noted that we have shifted some lines for clarity, but the slopes remain the same, which are the main concerns in our experiments. From Figures 5 and 6, we can conclude that the convergence rates in spatial and temporal directions coincide with error estimate in Theorem 1 are all of order 1.

5.   Conclusions

In this paper, we propose an efficient numerical method for the unilateral pricing problem of American better-of options with two underlying assets. By some traditional numeraire transformations, the far-field truncation technique and the known information about the free boundary, the original pricing model is approximated by a one-dimensional LCP on a bounded domain. In order to discretize the resulting LCP, the FDM and the FEM are applied in temporal direction and spatial direction, respectively. Based on the positive definite property of the discretized matrix, the discretized system is solved by the PDAS method. Our algorithm can obtain the price and the optimal exercise boundary of American better-of options simultaneously. Furthermore, we use some numerical experiments to verify the superiority of the proposed algorithm on error estimates and computational costs, respectively. In the future work, we shall apply our methods to different stochastic volatility models, regime switching models, and fractional models ^[32,33,34] to price American better-of options, as well as some inverse problems, such as the implied volatility of American options ^[35,36].

Acknowledgement

The work of K. Zhang was supported by the NSF of China under the grant No. 11871245, and by the Key Laboratory of Symbolic Computation and Knowledge Engineering of Ministry of Education, Jilin University (93K172018Z01). The work of H. Song was supported by the NSF of China under the grant No.11701210, the education department project of Jilin Province under the grant No. JJKH20211031KJ, the NSF of Jilin Province under the grants No.20190103029JH, 20200201269JC and the fundamental research funds for the Central Universities.

Conflict of interest

The authors declare there is no conflicts of interest.