Iteration changes discontinuity into smoothness (Ⅱ): oscillating case

Tianqi Luo; Xiaohua Liu; Tianqi Luo; Xiaohua Liu

doi:10.3934/math.2023441

AIMS Mathematics

2023, Volume 8, Issue 4: 8793-8810. doi: 10.3934/math.2023441

Previous Article Next Article

Research article

Iteration changes discontinuity into smoothness (Ⅱ): oscillating case

Tianqi Luo ,
Xiaohua Liu ^,

School of Mathematics and Physics, Leshan Normal University, Leshan 614000, China

Received: 15 April 2022 Revised: 26 June 2022 Accepted: 12 July 2022 Published: 08 February 2023
MSC : 37E05, 39B12

It has been shown that a self-mapping with exactly one removable or jumping discontinuity may have a $C^1$ smooth iterate of the second-order. However, some examples show that a self-mapping with exactly one oscillating discontinuity may also have a $C^1$ smooth iterate of the second-order, indicating that iteration can turn a self-mapping with exactly one oscillating discontinuity into a $C^1$ smooth one. In this paper, we study piecewise $C^1$ self-mappings on the open interval $(0, 1)$ having only one oscillating discontinuity. We give necessary and sufficient conditions for those self-mappings whose second-order iterates are $C^1$ smooth.

Keywords:

Citation: Tianqi Luo, Xiaohua Liu. Iteration changes discontinuity into smoothness (Ⅱ): oscillating case[J]. AIMS Mathematics, 2023, 8(4): 8793-8810. doi: 10.3934/math.2023441

Related Papers:

[1]	Min-Ku Lee, Jeong-Hoon Kim . Closed-form approximate solutions for stop-loss and Russian options with multiscale stochastic volatility. AIMS Mathematics, 2023, 8(10): 25164-25194. doi: 10.3934/math.20231284
[2]	Din Prathumwan, Thipsuda Khonwai, Narisara Phoochalong, Inthira Chaiya, Kamonchat Trachoo . An improved approximate method for solving two-dimensional time-fractional-order Black-Scholes model: a finite difference approach. AIMS Mathematics, 2024, 9(7): 17205-17233. doi: 10.3934/math.2024836
[3]	Boling Chen, Guohe Deng . Analytic approximations for European-style Asian spread options. AIMS Mathematics, 2024, 9(5): 11696-11717. doi: 10.3934/math.2024573
[4]	Zhongdi Cen, Jian Huang, Aimin Xu . A posteriori grid method for a time-fractional Black-Scholes equation. AIMS Mathematics, 2022, 7(12): 20962-20978. doi: 10.3934/math.20221148
[5]	Ashish Awasthi, Riyasudheen TK . An accurate solution for the generalized Black-Scholes equations governing option pricing. AIMS Mathematics, 2020, 5(3): 2226-2243. doi: 10.3934/math.2020147
[6]	Jafar Biazar, Fereshteh Goldoust . Multi-dimensional Legendre wavelets approach on the Black-Scholes and Heston Cox Ingersoll Ross equations. AIMS Mathematics, 2019, 4(4): 1046-1064. doi: 10.3934/math.2019.4.1046
[7]	Chao Yue, Chuanhe Shen . Fractal barrier option pricing under sub-mixed fractional Brownian motion with jump processes. AIMS Mathematics, 2024, 9(11): 31010-31029. doi: 10.3934/math.20241496
[8]	Yao Fu, Sisi Zhou, Xin Li, Feng Rao . Multi-assets Asian rainbow options pricing with stochastic interest rates obeying the Vasicek model. AIMS Mathematics, 2023, 8(5): 10685-10710. doi: 10.3934/math.2023542
[9]	Azadeh Ghanadian, Taher Lotfi . Approximate solution of nonlinear Black–Scholes equation via a fully discretized fourth-order method. AIMS Mathematics, 2020, 5(2): 879-893. doi: 10.3934/math.2020060
[10]	Kazem Nouri, Milad Fahimi, Leila Torkzadeh, Dumitru Baleanu . Numerical method for pricing discretely monitored double barrier option by orthogonal projection method. AIMS Mathematics, 2021, 6(6): 5750-5761. doi: 10.3934/math.2021339

Abstract

1. Introduction

The well-known Black-Scholes (B-S) equation ^[1] for multi-asset option pricing is the following d-dimensional partial differential equation (PDE):

$\begin{equation} \frac{\partial v({\bf S},t)}{\partial t} +\frac{1}{2}\sum\limits_{i,j = 1}^d\sigma_i \sigma_j \rho_{ij}S_i S_j\frac{\partial^2 v({\bf S},t)}{{\partial S_i}{\partial S_j} } +r \sum\limits_{i = 1}^d S_i \frac{\partial v({\bf S},t)}{\partial S_i} -r v({\bf S},t) = 0, \end{equation}$

(1.1)

for $({\bf S}, t) \in (0, \infty)^d\times [0, T)$ with the final condition $v({\bf S}, T) = v_T({\bf S})$ , where $v({\bf S}, t)$ is the value of the option in the multi-asset ${\bf S} = (S_1, S_2, \dots, S_d)$ at time $t$ . $T$ is the maturity time, $\sigma_i$ is the volatility of underlying asset $S_i$ , $\rho_{i j}$ is the correlation between $i$ -th and $j$ -th assets, and $r$ is the risk-free interest rate.

However most of the financial derivatives specially linear ones have analytical solutions, but they are not often easy to implement. Therefore the numerical methods would be suitable to find numerical approximations of them. Various numerical methods have been considered for one and multi-dimensional Black-Scholes equations. For instance, mesh free methods ^[9], high-order option pricing schemes ^[10,14], alternating direction implicit schemes ^[2,7] and etc. Furthermore, not only the numerical methods are so common to implement for the financial derivatives, but also they have been implemented for other variety problems and equations frequently ^[3,11,12,15].

One of the most common and simple numerical methods to implement for solving PDEs is FDMs. The main problem of FDMs for pricing multi-asset options with non-linear payoffs is significant numerical errors around their boundary points and also the strike price K. Several studies have been done to achieve an accurate numerical solution of multi-asset options. Jeong et al. ^[8] proposed a useful remedy to reduce numerical error around ${\bf S}_{max}$ by using the standard Monte Carlo simulation at this boundary point. Here we improved their methodology by using the antithetic variate method in Monte Carlo simulation ^[4], and furthermore a space transformation ^[13,16] that causes more nodal points around the strike price to get more robust solution.

2. FDM for the two-asset B-S equation with combination of the GS and the MC simulation

Two-dimensional B-S equation with two assets $x = S_1$ and $y = S_2$ is the following PDE:

$\begin{equation} \frac{\partial u}{\partial \tau} = \frac{1}{2}\sigma^2_x x^2\frac{\partial^2 u}{\partial x^2} +\frac{1}{2}\sigma_y^2 y^2\frac{\partial^2 u}{\partial y^2} +\sigma_x\sigma_y \rho\frac{\partial^2 u}{{\partial x}{\partial y} } +r x \frac{\partial u}{\partial x} +r y \frac{\partial u}{\partial y} -r u , \end{equation}$

(2.1)

for $(x, y, \tau) \in \Omega\times [0, T)$ with the initial condition $u(x, y, 0) = u_0(x, y)$ which is the Eq (1.1) with $d = 2$ and a variable changing $\tau = T-t$ in the truncated domain $\Omega = (0, x_{max})\times(0, y_{max})$ . Suppose $\mathcal{L}_{BS}$ is the operator of the right hand side of Eq (2.1). Now we discretize the time interval with a uniform time step $\Delta \tau = T/N_{\tau}$ and the following grid stretching transformation that makes non-uniform space discretization of the interval $[x_{min}, x_{max}]$ with more nodal points around $K$ :

$\begin{equation} x = K\left(1+\frac{\sinh(c_2w+c_1(1-w))}{15}\right) , \ \ {\rm for } \ \ w\in[0,1], \end{equation}$

(2.2)

where $c_1 = {\rm{arcsinh}}(15(x_{min}-K)/K)$ and $c_2 = {\rm{arcsinh}}(15(x_{max}-K)/K)$ . Therefore $x(w)\in[x_{min}, x_{max}]$ for any $w\in[0, 1]$ . The grid stretching transformation by locating more nodal points around the strike price $K$ , leads to reduce the numerical errors around $K$ for the options with non-smooth payoffs which can be applied to even higher order finite difference methods in both single and multi-asset Black-Scholes equations. At first for the $x$ -direction, we discretize the interval $[0, 1]$ uniformly to $N_x$ sub-intervals, then we use the transformation (2.2) in $[0, x_{max}]$ and define $h_i^x = x_{i+1}-x_i$ for $i = 0, \dots, N_x-1$ . we discretize the interval $[0, y_{max}]$ to $N_y$ sub-intervals similarly. Let for $i = 0, \dots, N_x$ , $j = 0, \dots, N_y$ and $n = 0, \dots, N_{\tau}$ define $u^n_{ij} \approx u(x_i, y_j, \tau_n)$ and use the operator splitting method for solving Eq (2.1) with the following two discrete fractional time step equations:

$\begin{equation} \frac{u_{ij}^{n+\frac{1}{2}}-u_{ij}^n}{\Delta \tau} = \mathcal{L}_{BS}^x\left(\alpha u_{ij}^{n+\frac{1}{2}} + \beta u_{ij}^{n}\right), \end{equation}$

(2.3)

$\begin{equation} \frac{u_{ij}^{n+1}-u_{ij}^{n+\frac{1}{2}}}{\Delta \tau} = \mathcal{L}_{BS}^y\left(\alpha u_{ij}^{n+1} + \beta u_{ij}^{n+\frac{1}{2}}\right), \end{equation}$

(2.4)

where $\alpha \geq 0$ , $\beta \geq 0$ and $\alpha+\beta = 1$ . We can choose different sets of values of $\alpha$ and $\beta$ to achieve us implicit or explicit schemes in both (2.3) and (2.4). These schemes have first-order accuracy in general ^[5]. Here we take the implicit scheme by choosing $\alpha = 1$ and $\beta = 0$ for both Eqs (2.3) and (2.4) which leads to a consistent and stable scheme and thereby the scheme is convergent with the first order of accuracy in time. Furthermore, with the operators $\left(\mathcal{L}_{BS}^x u\right)_{ij}^{n+\frac{1}{2}}$ and $\left(\mathcal{L}_{BS}^y u\right)_{ij}^{n+1}$ as follows, we obtain the second order of accuracy in space as the central finite difference approximations have been used for the space derivatives.

$\begin{equation} \left(\mathcal{L}_{BS}^x u\right)_{ij}^{n+\frac{1}{2}} = \frac{(\sigma_x x_i)^2}{2}D_{xx}u_{ij}^{n+\frac{1}{2} }+ \frac{\sigma_x \sigma_y \rho x_i y_j }{2}D_{xy} u_{ij}^n + r x_i D_x u_{ij}^{n+\frac{1}{2} } -\frac{r}{2} u_{ij}^{n+\frac{1}{2} }, \end{equation}$

(2.5)

$\begin{equation} \left(\mathcal{L}_{BS}^y u\right)_{ij}^{n+1} = \frac{(\sigma_y y_j)^2}{2}D_{yy}u_{ij}^{n+1} + \frac{\sigma_x \sigma_y \rho x_i y_j }{2}D_{xy} u_{ij}^{n+\frac{1}{2}}+ r y_j D_y u_{ij}^{n+1 } -\frac{r}{2} u_{ij}^{n+1}, \end{equation}$

(2.6)

where the derivatives approximations are defined as

$\begin{equation} D_{xx} u_{ij} \approx \frac{2}{h^x_{i-1}\;(h^x_{i-1}\;+h^x_{i})} u_{i-1,j} - \frac{2}{h^x_{i-1}\;h^x_{i}} u_{ij} + \frac{2}{h^x_{i}(h^x_{i-1}\;+h^x_{i})} u_{i+1,j} , \end{equation}$

(2.7)

$\begin{equation} D_{xy} u_{ij} \approx \frac{u_{i+1,j+1}\; - u_{i-1,j+1}\; - u_{i+1,j-1}\; + u_{i-1,j-1}}{h^x_{i} h^y_{j}+h^x_{i-1}\;h^y_{j} + h^x_{i} h^y_{j-1}\;+ h^x_{i-1}\;h^y_{j-1}} , \end{equation}$

(2.8)

$\begin{equation} D_{x} u_{ij} \approx -\frac{h_i^x}{h^x_{i-1}\;(h^x_{i-1}\;+h^x_{i})} u_{i-1,j} + \frac{h_i^x - h_{i-1}^x}{h^x_{i-1}\;h^x_{i}} u_{ij} + \frac{h_{i-1}^x}{h^x_{i}(h^x_{i-1}\;+h^x_{i})} u_{i+1,j}. \end{equation}$

(2.9)

By substituting above derivatives in Eq (2.3), the following equations will be achieved:

$\begin{equation} \alpha_i u^{n+\frac{1}{2}}_{i-1,j} = \beta_i u^{n+\frac{1}{2}}_{i,j} + \gamma_i u^{n+\frac{1}{2}}_{i+1,j} = f_{ij} , \ \ \ i = 1,\dots , N_x , \end{equation}$

(2.10)

where $\alpha_i$ , $\beta_i$ and $\gamma_i$ are the lower, main and upper diagonal of the tridiagonal matrix of the above system of equations respectively and $f_{ij}$ is the constant vector of the system which are defined as:

$\begin{equation} \alpha_i = -\frac{(\sigma_x x_i)^2}{h_{i-1}^x\;(h_{i-1}^x\;+h_{i}^x)} + \frac{r x_i h_i^x}{h_{i-1}^x\;(h_{i-1}^x\;+h_{i}^x)} , \end{equation}$

(2.11)

$\begin{equation} \beta_i = \frac{1}{\Delta \tau} + \frac{(\sigma_x x_i)^2}{h_{i-1}^x\;h_{i}^x} - \frac{r x_i (h_i^x\;- h_{i-1}^x\;)}{h_{i-1}^x\;h_{i}^x} + \frac{r}{2} , \end{equation}$

(2.12)

$\begin{equation} \gamma_i = - \frac{(\sigma_x x_i)^2}{h_{i}^x(h_{i-1}^x\;+h_{i}^x)} - \frac{r x_i h_{i-1}^x}{h_{i}^x(h_{i-1}^x\;+h_{i}^x)} , \end{equation}$

(2.13)

$\begin{equation} f_{ij}^n = \frac{1}{2}\sigma_x \sigma_y \rho x_i y_j D_{xy} u_{ij}^n - \frac{1}{\Delta \tau} u_{ij}^n, \end{equation}$

(2.14)

for fixed $j$ and the linear boundary condition at the boundaries. For instance $u_{0j} = 2u_{1j}-u_{2j}$ , for $j = 1, \dots, N_y$ and similarly for $u_{N_x, j}$ , $u_{i0}$ and $u_{i, N_y}$ or at the end point of the domain, we apply the Monte Carlo simulation with the antithetic variate method (AMC) or the standard Monte Carlo simulation (MC). The solution of the above system of equations will be the vector $u_{1:N_x, j}^{\frac{1}{2}}$ . Similarly by substituting the derivatives approximations with respect to $y$ in (2.4), we have another system of equations to solve for finding the value of the two-asset option.

Note that with the AMC and the MC simulations we obtain the option value $V_{AMC}$ and $V_{MC}$ respectively at time $t = 0$ and $(x, y) = (x_{max}, y_{max})$ as follow:

$\begin{equation} V_{AMC} = \frac{e^{-rT}}{M}\sum\limits_{m = 1}^M\left[ \frac{1}{2} \left( \Lambda({\bf S}_{max}) + \Lambda(\tilde{ {\bf S}}_{max} ) \right) \right] , \end{equation}$

(2.15)

$\begin{equation} V_{MC} = \frac{e^{-rT}}{M}\sum\limits_{m = 1}^M \left( \Lambda({\bf S}_{max}) \right) , \end{equation}$

(2.16)

where

$\begin{equation} {\bf S}_{max} = \left[ x_{max}\times e^{\left((r-\frac{\Phi^2}{2})T+\Sigma(1,:).{\bf z}.\sqrt{T}\right)} , y_{max}\times e^{\left((r-\frac{\Phi^2}{2})T+\Sigma(2,:).{\bf z}.\sqrt{T}\right) }\right] , \end{equation}$

(2.17)

$\begin{equation} \tilde{{\bf S}}_{max} = \left[ x_{max}\times e^{\left((r-\frac{\Phi^2}{2})T-\Sigma(1,:).{\bf z}.\sqrt{T}\right)} , y_{max}\times e^{\left((r-\frac{\Phi^2}{2})T-\Sigma(2,:).{\bf z}.\sqrt{T}\right) }\right], \end{equation}$

(2.18)

and $M$ is the number of replications in the antithetic and the standard Monte Carlo simulation which in this work we choose $M = 10^6$ , $\Lambda$ is a payoff function, ${\bf z}$ , $\Sigma$ and $\Phi$ are the random vector from a standard normal distribution, the volatility matrix and a vector with the size of multi-asset respectively. For instance, for the case of two-asset $\Sigma$ and $\Phi$ are as follows:

$\begin{equation} \Sigma = \left[ {\begin{array}{cc} \sigma_{11} & \sigma_{12} \\ \sigma_{21} & \sigma_{22} \end{array} } \right], \end{equation}$

(2.19)

$\begin{equation} \Phi = \sqrt{diag(\Sigma \times \Sigma')} = \left[ \begin{array}{c} \sqrt{\sigma_{11}\times \sigma_{12}} \\ \sqrt{\sigma_{21}\times \sigma_{22}} \end{array} \right] , \end{equation}$

(2.20)

where the elements of the volatility matrix will be obtained by solving the following system of equations:

$\begin{eqnarray} \sigma_x = \sqrt{\sigma^2_{11}+\sigma^2_{12}}, \end{eqnarray}$

(2.21)

$\begin{eqnarray} \sigma_y = \sqrt{\sigma^2_{21}+\sigma^2_{22}}, \end{eqnarray}$

(2.22)

$\begin{eqnarray} \rho = \frac{\sigma_{11}\sigma_{21}+\sigma_{12}\sigma_{22}}{\sqrt{(\sigma^2_{11}+\sigma^2_{12})(\sigma^2_{21}+\sigma^2_{22})}}. \end{eqnarray}$

(2.23)

Then with the exponential interpolation between $V$ (which is $V_{AMC}$ or $V_{MC}$ ) and the payoff value $u_{N_x, N_y}^0$ , we obtain the option price at the end point of the domain for every time step:

$\begin{equation} u_{N_x,N_y}^n = u_{N_x,N_y}^0 \left( \frac{V}{u_{N_x,N_y}^0} \right)^{\frac{n\Delta \tau}{T}} , \quad \ \ n = 0,\dots, N_{\tau}. \end{equation}$

(2.24)

3. Numerical results

Here we solve numerically the B-S equation in one and two-dimensional space for different options with these parameters: $T = 1$ , $K_1 = K_2 = 100$ , $r = 0.03$ , $\sigma_x = \sigma_y = 0.3$ , $\rho = 0.5$ .

3.1. Power option

A power call option has a payoff $max(x^p-K, 0)$ where $p\in \mathcal{R}^+$ is a power. The closed-form solution of the power option is given by ^[18]:

$\begin{equation} u(x,\tau) = x^pe^{(p-1)(r+p\sigma^2/2)\tau} N(d_1)-Ke^{-r\tau}N(d_2) \end{equation}$

(3.1)

where

$\begin{eqnarray} d_1 & = \frac{\ln \left( \frac{x}{\sqrt[p]{K}} \right)+ \left( r+ (p-0.5)\sigma^2\right) \tau }{\sigma \sqrt{\tau}}\, , \qquad d_2 = d_1-p\sigma \sqrt{\tau}. \end{eqnarray}$

(3.2)

Note that the power option with $p = 1$ is the call option. shows the comparison of the numerical errors of one-asset power call option with different schemes: The equidistance discretization in space with the linear boundary condition (Equi&L), the grid stretching discretization around the strike price with linear boundary condition (GS& $L$ ), the equidistance discretization in space with Monte Carlo simulation at $S_{max}$ (Equi&MC) and the grid stretching discretization with the antithetic Monte Carlo simulation at $S_{max}$ (GS&AMC).

Figure 1. One-asset power call option with

$p = 2$ ,

$S_{max} = \sqrt{250}$ ,

$N_x = 250$ and

$\Delta \tau = \frac{1}{720}$ , exact option price (top-left), payoff (top-right), Equi&L error (middle-left), GS&L error (middle-right), Equi&MC error (bottom-left) and GS&AMC error (bottom-right).

DownLoad: Full-Size Img PowerPoint

shows the RMSE of the power option with $N_x = 250$ and $\Delta \tau = \frac{1}{720}$ for $p = 1$ , $p = 2$ and $p = 3$ in the whole domain $\Omega = (0, S_{max})$ and the interesting domain around the strike price $\Omega_e = [0.7\sqrt[p]{K}, 1.3\sqrt[p]{K}]$ for the different schemes and different $S_{max}$ . Keep in mind that in the grid stretching transformation (2.2) we must substitute $\sqrt[p]{K}$ instead of $K$ . We can see the MC and AMC simulation at the end point of the domain reduce RMSE significantly for bigger $p$ and none of the GS, MC and AMC improve the accuracy of the numerical approximation for the call option ( $p = 1$ ).

Table 1. RMSE for the power call option.

	$S_{max}$	150	200	250	300	350	400	450
	Equi&L	0.3452	0.0271	0.0024	0.0009	0.0008	0.0008	0.0008
p=1	GS&L	0.2737	0.0165	0.0022	0.0017	0.0016	0.0016	0.0015
$\Omega$	Equi&MC	0.0457	0.0172	0.0072	0.0511	0.0148	0.0038	0.0419
	GS&AMC	0.0457	0.0033	0.0032	0.0023	0.0023	0.0017	0.0024
	Equi&L	0.1505	0.0015	0.0015	0.0016	0.0017	0.0018	0.0019
p=1	GS&L	0.1122	0.0027	0.0024	0.0022	0.0021	0.0021	0.0020
$\Omega_K$	Equi&MC	0.0563	0.0016	0.0015	0.0016	0.0017	0.0018	0.0019
	GS&AMC	0.0472	0.0023	0.0024	0.0022	0.0021	0.0021	0.0020
	$S_{max}$	$\sqrt{150}$	$\sqrt{200}$	$\sqrt{250}$	$\sqrt{300}$	$\sqrt{350}$	$\sqrt{400}$	$\sqrt{450}$
	Equi&L	6.8424	6.8064	7.7089	8.9340	10.2820	11.6811	13.1029
p=2	GS&L	8.4653	6.1009	5.8106	6.0578	6.4852	6.9842	7.5149
$\Omega$	Equi&MC	0.2272	0.0685	0.0319	0.0205	0.0457	0.0926	0.0863
	GS&AMC	0.3311	0.0253	0.0244	0.0345	0.0266	0.0364	0.0348
	Equi&L	10.4749	5.5909	2.9152	1.7040	1.0414	0.6631	0.4350
p=2	GS&L	9.6261	4.6061	2.3610	1.3371	0.8022	0.5133	0.3308
$\Omega_K$	Equi&MC	0.3478	0.0661	0.0193	0.0170	0.0115	0.0161	0.0106
	GS&AMC	0.3765	0.0245	0.0165	0.0203	0.0146	0.0117	0.0081
	$S_{max}$	$\sqrt[3]{150}$	$\sqrt[3]{200}$	$\sqrt[3]{250}$	$\sqrt[3]{300}$	$\sqrt[3]{350}$	$\sqrt[3]{400}$	$\sqrt[3]{450}$
	Equi&L	19.3431	22.0844	25.8076	29.9848	34.3954	38.9355	43.5521
p=3	GS&L	28.5178	25.4178	25.0462	25.8940	27.3175	29.0375	30.9201
$\Omega$	Equi&MC	1.2650	0.6127	0.4785	0.3064	0.3175	0.2029	0.1936
	GS&AMC	1.9696	0.9038	0.5715	0.3920	0.3462	0.2689	0.2511
	Equi&L	31.1024	33.1054	28.4051	21.9288	17.3308	14.1742	11.8378
p=3	GS&L	32.8550	28.7798	23.0077	17.4396	13.9633	11.2298	9.3862
$\Omega_K$	Equi&MC	2.0303	0.9183	0.6883	0.4127	0.3427	0.2148	0.1607
	GS&AMC	2.2680	1.0233	0.6413	0.4100	0.3167	0.2264	0.1849

| Show Table

DownLoad: CSV

3.2. European call option on the maximum of two assets

Now we consider a European call option on the maximum of two assets with the payoff $u(x, y, 0) = max\{max\{x, y\}-K, 0\}$ . Then the following exact analytical solutions of the option ^[17] will be compared to the numerical solutions.

$\begin{align} u(x,y,\tau) = & x M\left(d_1,d;\rho_1\right) + y M\left(d_2,-d+\sigma \sqrt {\tau};\rho_2\right) \\ &- K e^{-r\tau}\left[ 1-M\left(-d_1+\sigma_1\sqrt{\tau}, -d_2+\sigma_2\sqrt{\tau}; \rho \right) \right] , \end{align}$

(3.3)

where $M$ is the cumulative bivariate normal distribution function and defined as

$\begin{equation} M(a,b;\rho) = \frac{1}{2\pi\sqrt{1-\rho^2}}\int_{-\infty}^b\int_{-\infty}^a \exp\left( - \frac{x^2-2\rho x y +y^2}{2(1-\rho^2)}\right)\, dx dy\, , \end{equation}$

(3.4)

and other parameters are as follows:

$\begin{eqnarray} d_1 & = \frac{\ln \left( \frac{x}{K} \right)+ \left( r+ \frac{\sigma_x^2}{2}\right) T }{\sigma_x\sqrt{T}}\, , d_2 = \frac{\ln \left( \frac{y}{K} \right)+ \left( r+ \frac{\sigma_y^2}{2}\right) T }{\sigma_y\sqrt{T}}\, , d = \frac{\ln \left( \frac{x}{y} \right)+ \frac{\sigma^2}{2} T }{\sigma\sqrt{T}}\, , \end{eqnarray}$

(3.5)

$\begin{eqnarray} \sigma & = \sqrt{\sigma_x^2 + \sigma_y^2 -2\rho \sigma_x\sigma_y }\, , \ \ \ \ \rho_1 = (\sigma_x - \rho \sigma_y )/\sigma \, , \ \ \ \ \rho_2 = (\sigma_y - \rho \sigma_x )/\sigma. \end{eqnarray}$

(3.6)

shows the comparison of the numerical solutions of the call option on the maximum of two assets with the different schemes and $(x_{max}, y_{max}) = (250,250)$ .

Figure 2. Call option on the maximum of two assets with

$N_x = Ny = 250$ and

$\Delta \tau = \frac{1}{720}$ exact option price (top-left), payoff (top-right), Equi&L error (middle-left), GS&L error (middle-right), Equi&MC error (bottom-left) and GS&AMC error (bottom-right).

DownLoad: Full-Size Img PowerPoint

Now in we compare the RMSE of these four schemes in the whole region $\Omega = (0, x_{max})\times (0, y_{max})$ and in the most interesting region $\Omega_K = (0.7K, 1.3K)\times(0.7K, 1.3K)$ with different domain sizes $(x_{max}, y_{max}) = (L, L)$ and $N_x = N_y = L$ at $t = 0$ .

Table 2. RMSE of call on the maximum of two assets.

	L	150	200	250	300
	Equi&L	12.1428	16.9562	22.3907	28.4710
	GS&L	9.1971	5.7377	5.4052	5.5654
$\Omega$	Equi&MC	2.6188	3.0868	3.8931	4.6910
	GS&AMC	2.2287	1.5009	1.7865	2.0569
	Equi&L	8.5065	0.7793	0.0690	0.0084
	GS&L	6.1609	0.4476	0.0354	0.0051
$\Omega_K$	Equi&MC	2.1547	0.2287	0.0225	0.0048
	GS&AMC	1.9339	0.1858	0.0178	0.0040

| Show Table

DownLoad: CSV

3.3. European call on the minimum of two assets

A European call on the minimum of two assets has the payoff $u(x, y, 0) = max\{min\{x, y\}-K, 0\}$ and the closed form solution as follows ^[17]

$\begin{align} u(x,y,\tau) = & x M\left(d_1,-d;-\rho_1\right) + y M\left(d_2,d-\sigma \sqrt {\tau};-\rho_2\right) \\ &- K e^{-r\tau}M\left(d_1-\sigma_1\sqrt{\tau}, d_2-\sigma_2\sqrt{\tau}; \rho \right). \end{align}$

(3.7)

where its parameters are the same as the European call option on the maximum of two assets. In Figure 3 and Table 3 we compare the four above mentioned schemes for this two-asset option.

Figure 3. Call option on the minimum of two assets with

$(x_{max}, y_{max}) = (250,250)$ ,

$N_x = Ny = 250$ and

$\Delta \tau = \frac{1}{720}$ exact option price (top-left), payoff (top-right), Equi&L error (middle-left), GS&L error (middle-right), Equi&MC error (bottom-left) and GS&AMC error (bottom-right).

DownLoad: Full-Size Img PowerPoint

Table 3. RMSE of call on the minimum of two assets.

	L	150	200	250	300
	Equi&L	11.4121	16.5529	21.9750	28.0137
	GS&L	8.7458	5.6337	5.1434	5.1969
$\Omega$	Equi&MC	0.9972	1.4289	1.8595	3.7692
	GS&AMC	0.8470	0.7205	0.6078	0.5893
	Equi&L	8.2864	0.7763	0.0666	0.0062
	GS&L	5.9848	0.4443	0.0330	0.0030
$\Omega_K$	Equi&MC	1.0381	0.1382	0.0128	0.0029
	GS&AMC	0.6490	0.0951	0.0087	0.0016

| Show Table

DownLoad: CSV

3.4. Two-asset correlation call option

Here we consider a two-asset correlation call option wich pays off:

$\begin{equation} u(x,y,0) = \left\{ \begin{array}{cc} y-K_2 \ \ , & if \ \ x > K_1 \\ 0 \ \ , & otherwise \end{array} \right. \end{equation}$

(3.8)

This option has been priced by ^[18]:

$\begin{equation} u(x,y,\tau) = y M\left(d_y+\sigma_y\sqrt{\tau},d_x+\rho \sigma_y \sqrt {\tau};\rho\right) - Ke^{-r\tau}M\left(d_y, d_x; \rho \right), \end{equation}$

(3.9)

where $\rho$ is the correlation coefficient between the returns on the two assets and

$\begin{eqnarray} d_x & = \frac{\ln \left( \frac{x}{K} \right)+ \left( r- \frac{\sigma_x^2}{2}\right) T }{\sigma_x\sqrt{T}}\, , \ \ \ \ \ d_y = \frac{\ln \left( \frac{y}{K} \right)+ \left( r- \frac{\sigma_y^2}{2}\right) T }{\sigma_y\sqrt{T}}\, . \end{eqnarray}$

(3.10)

Now we solve numerically this option. Figure 4 shows the comparison of the numerical solutions of two-asset correlation call option with the different schemes.

Figure 4. Two-asset correlation call option with

$(x_{max}, y_{max}) = (250,250)$ ,

$N_x = Ny = 250$ and

$\Delta \tau = \frac{1}{720}$ , exact option price (top-left), payoff (top-right), Equi&L error (middle-left), GS&L error (middle-right), Equi&MC error (bottom-left) and GS&AMC error (bottom-right).

DownLoad: Full-Size Img PowerPoint

shows the RMSE of the equidistance discretization (Equi) and the grid stretching discretization around the strike price with Monte Carlo simulation at the last point of the boundary $(x_{max}, y_{max}) = (L, L)$ (GS&MC)with different $L$ in the whole region $\Omega = [0, L]\times[0, L]$ and the most interesting region $\Omega_K = [0.7K, 1.3K]\times[0.7K, 1.3K]$ respectively.

Table 4. RMSE for the correlation call option.

	L	150	200	250	300
	Equi&L	0.8040	0.5688	0.6639	0.7581
	GS&L	0.7340	0.6739	0.9008	0.9729
$\Omega$	Equi&MC	0.2578	0.3740	0.5793	0.8419
	GS&AMC	0.1999	0.5024	0.8245	0.8788
	Equi&L	0.2386	0.1022	0.1023	0.1023
	GS&L	0.1803	0.0213	0.0181	0.0108
$\Omega_K$	Equi&MC	0.2582	0.1255	0.1047	0.1026
	GS&AMC	0.2355	0.0474	0.0171	0.0110

| Show Table

DownLoad: CSV

3.5. Two-asset butterfly option

A butterfly option can be regarded as a combination of two long calls with strikes $K_1$ and $K_2$ and two short calls both with strike $K = (K_1 + K_2)/2$ . Thus, for given values $K_1, K_2 > 0$ with $K_1 < K_2$ , the payoff of a two-asset European butterfly option is:

$\begin{eqnarray} u(x,y,0) = &max\{max\{x,y\}-K_1,0\}+max\{max\{x,y\}-K_2,0\}\\ &-2 max\{max\{x,y\}-K,0\}, \end{eqnarray}$

(3.11)

and has the following analytical solution:

$\begin{align*} \label{butterfly} u(x,y,\tau)& = x\left[ M\left(d_1^{(1)},d;\rho_1\right)-2M\left(d_1,d;\rho_1\right)+M\left(d_1^{(2)},d;\rho_1\right)\right] \\ &+ y\left[M\left(d_2^{(1)},-d+\sigma \sqrt {\tau};\rho_2\right)-2M\left(d_2,-d+\sigma \sqrt {\tau};\rho_2\right)+M\left(d_2^{(2)},-d+\sigma \sqrt {\tau};\rho_2\right)\right] \nonumber \\ &-K_1 e^{-r\tau}\left[ 1-M\left(-d_1^{(1)}+\sigma_x\sqrt{\tau}, -d_2^{(1)}+\sigma_y\sqrt{\tau}; \rho \right) \right] \\ &+2K e^{-r\tau}\left[ 1-M\left(-d_1+\sigma_x\sqrt{\tau}, -d_2+\sigma_y\sqrt{\tau}; \rho \right) \right] \\ &-K_2 e^{-r\tau}\left[ 1-M\left(-d_1^{(2)}+\sigma_x\sqrt{\tau}, -d_2^{(2)}+\sigma_y\sqrt{\tau}; \rho \right) \right], \nonumber \end{align*}$

where the $M(a, b;\rho)$ is the cumulative bivariate normal distribution function (3.4) and other parameters are defined in (3.5) and (3.6). where $d_1^{(i)}$ and $d_2^{(i)}$ will be obtained by substituting $K_i$ in 3.5 for $i = 1, 2$ .

Here we consider a two-asset butterfly option with $K_1 = 50$ , $K_2 = 150$ . shows the comparison of the numerical solutions of the two-asset butterfly option with the equidistance and the grid stretching discretization schemes. Since $u^0_{N_x, N_y} = 0$ we can not use the Monte Carlo simulation at the end point of the domain in (2.24).

Figure 5. Two-asset butterfly option with

$K_1 = 50, K_2 = 150$ ,

$(x_{max}, y_{max}) = (250,250)$ ,

$N_x = Ny = 250$ and

$\Delta \tau = \frac{1}{720}$ , exact option price (top-left), payoff (top-right), Equi error (bottom-left), GS error (bottom-right).

DownLoad: Full-Size Img PowerPoint

Now in we compare the RMSE of the equidistance discretization (Equi) and the grid stretching discretization (GS) for the space in the whole region $\Omega = (x_{max}, y_{max})$ and the most interesting region $\Omega_K = [0.7K_1, 1.3K_2]\times[0.7K_1, 1.3K_2]$ with different domain sizes $(x_{max}, y_{max}) = (L, L)$ and $N_x = N_y = L$ at $t = 0$ . We can see in this option the GS scheme does not have any superiority over Equi for the bigger domain.

Table 5. RMSE for the two-asset butterfly option.

	L	200	250	300	350
	Equi	1.0136	0.2986	0.2065	0.1753
$\Omega$	GS	0.5842	0.2733	0.2258	0.1946
	Equi	0.8978	0.0174	0.0043	0.0042
$\Omega_K$	GS	0.5099	0.0098	0.0052	0.0049

| Show Table

DownLoad: CSV

3.6. Two-asset cash-or-nothing call option

The payoff function of a two-asset cash-or-nothing option is given by

$\begin{equation} u(x,y,0) = \left\{ \begin{array}{cc} c , & if \ \ x\geq K_1 , y\geq K_2 \\ 0 , & otherwise, \end{array} \right. \end{equation}$

(3.12)

which $c$ is a constant. The closed form solution for the two-asset cash-or-nothing option is as follows

$\begin{equation} u(x,y,\tau) = ce^{-r\tau}M(d_x,d_y;\rho), \end{equation}$

(3.13)

where $M(a, b;\rho)$ is the cumulative bivariate normal distribution function (3.4) and $d_x$ and $d_y$ are the same as (3.10). Now we consider a two-asset cash-or-nothing call option with $c = 100$ , then the exact analytical solutions of the option ^[6] will be compared to the numerical solutions. shows the comparison of the numerical solutions of the option with the different four schemes and also demonstrates the RMSE of the four schemes with different domain sizes $(x_{max}, y_{max}) = (L, L)$ , $N_x = N_y = L$ and $\Delta \tau = \frac{1}{720}$ in the whole region $\Omega$ and the in interesting region $\Omega_K$ . Since the payoff of this option is constant at the end point of the domain, we can see that the MC and AMC simulations do not have any superiority for the bigger domain.

Figure 6. Two-asset cash-or-nothing call option with

$(x_{max}, y_{max}) = (250,250)$ ,

$N_x = Ny = 250$ and

$\Delta \tau = \frac{1}{720}$ , exact option price (top-left), payoff (top-right), Equi&L error (middle-left), GS&L error (middle-right), Equi&MC error (bottom-left) and GS&AMC error (bottom-right).

DownLoad: Full-Size Img PowerPoint

Table 6. RMSE for the cash-or-nothing call option.

	L	150	200	250	300
	Equi&L	3.6561	1.1174	0.9870	0.9458
	GS&L	3.0292	0.6828	0.5302	0.3400
$\Omega$	Equi&MC	1.3850	1.0302	0.9847	0.9462
	GS&AMC	1.1707	0.6454	0.5294	0.3400
	Equi&L	0.0429	0.0218	0.0217	0.0217
	GS&L	0.0237	0.0051	0.0042	0.0024
$\Omega_K$	Equi&MC	0.0250	0.0217	0.0217	0.0217
	GS&AMC	0.0047	0.0050	0.0042	0.0024

| Show Table

DownLoad: CSV

4. Conclusions

We applied the FDM for solving numerically one and two-dimensional B-S equations with different space discretization and boundary conditions. Also, RMSE has been applied for estimating error of the method as it is common to use and consider the numerical errors for all nodal points of the region. Numerical experiments illustrated that by applying the antithetic Monte Carlo simulation for the end point of the space domain and discretization with grid stretching to put more nodal points around the non-smooth point of the payoff function, we obtained more accurate numerical solutions in the considered exotic options with non-constant payoff at the end point of the space domain. These remedies can be considered and implemented for other options, higher dimensional models, nonlinear B-S equations and also other financial derivatives.

Acknowledgments

The authors would like to thank the reviewers for their helpful suggestions and comments, which improved the quality of this article.

Conflict of interest

The authors declare that they have no conflicts of interest.

References

[1]	N. H. Abel, Oeuvres complétes, Vol. Ⅱ, Christiana, 1881, 36–39.
[2]	J. M. Dubbey, The mathematical work of Charles Babbage, Cambridge: Cambridge University Press, 1978. https://doi.org/10.1017/CBO9780511622397
[3]	X. Liu, L. Liu, W. Zhang, Discontinuous function with continuous second iterate, Aequat. Math., 88 (2014), 243–266. https://doi.org/10.1007/s00010-013-0220-z doi: 10.1007/s00010-013-0220-z
[4]	X. Liu, L. Liu, W. Zhang, Smoothness repaired by iteration, Aequat. Math., 89 (2015), 829–848. https://doi.org/10.1007/s00010-014-0277-3 doi: 10.1007/s00010-014-0277-3
[5]	X. Liu, S. Luo, Continuity of functions with finitely many discontinuities of the same type repaired by iteration, (Chinese), Acta Math. Sinica, 65 (2022), 123–146. https://doi.org/10.12386/A20220010 doi: 10.12386/A20220010
[6]	T. Luo, X. Liu, Iteration changes discontinuity into smoothness (Ⅰ): removable and jumping cases, AIMS Mathematics, 8 (2023), 8772–8792. https://doi.org/10.3934/math.2023440 doi: 10.3934/math.2023440
[7]	E. Schröder, Über iterate funktionen, Math. Ann., 3 (1871), 296–322.

This article has been cited by:

Y. Ma, C.Z. Shi, Y.C. Hon, Generalized Finite Integration Method with Laplace transform for European option pricing under Black–Scholes and Heston models, 2024, 164, 09557997, 105751, 10.1016/j.enganabound.2024.105751

Reader Comments

Your name:*

Email:*
© 2023 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)