The instability of periodic solutions for a population model with cross-diffusion

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):

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:

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$:

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:

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.

where the derivatives approximations are defined as

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

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:

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:

where

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:

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

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:

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]:

where

Note that the power option with $p = 1$ is the call option. Figure 1 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).

Table 1 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$).

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.

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

and other parameters are as follows:

Figure 2 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)$.

Now in Table 2 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$.

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]

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.

3.4.   Two-asset correlation call option

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

This option has been priced by ^[18]:

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

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.

Table 4 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.

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:

and has the following analytical solution:

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$. Figure 5 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).

Now in Table 5 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.

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

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

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

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. Figure 6 shows the comparison of the numerical solutions of the option with the different four schemes and also Table 6 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.

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.