Pricing foreign exchange options under stochastic volatility and interest rates using an RBF--FD method

TL;DR

This paper introduces a localized RBF-FD numerical method for efficiently pricing FX options under a complex model with stochastic interest rates and volatility, overcoming computational challenges of traditional methods.

Contribution

The paper develops a novel RBF-FD approach tailored for a complex FX option pricing PDE with non-affine terms, improving efficiency and accuracy.

Findings

The RBF-FD method achieves high accuracy in FX option pricing.

The approach reduces computational cost compared to traditional methods.

Numerical results demonstrate the method's effectiveness in computing Greeks.

Abstract

This paper proposes a numerical method for pricing foreign exchange (FX) options in a model which deals with stochastic interest rates and stochastic volatility of the FX rate. The model considers four stochastic drivers, each represented by an It\^{o}'s diffusion with time--dependent drift, and with a full matrix of correlations. It is known that prices of FX options in this model can be found by solving an associated backward partial differential equation (PDE). However, it contains non--affine terms, which makes its difficult to solve it analytically. Also, a standard approach of solving it numerically by using traditional finite--difference (FD) or finite elements (FE) methods suffers from the high computational burden. Therefore, in this paper a flavor of a localized radial basis functions (RBFs) method, RBF--FD, is developed which allows for a good accuracy at a relatively low computational cost. Results of numerical simulations are presented which demonstrate efficiency of such an approach in terms of both performance and accuracy for pricing FX options and computation of the associated Greeks.

Tables

Table 1. Table 1 : Parameters of the model used for the numerical experiments 1 and 2 .

$T,\text{yr}$	$E,\$$	$\gamma$	$\kappa$	$\overline{v}$	${\eta }_d$	${\eta }_f$	$r_{d,0}$	$r_{f,0}$
1,2	100	0.3	0.5	0.1	0.007	0.012	0.1	0.1
$v_0$	${\lambda }_d$	${\lambda }_f$	${\varrho }_1$	${\varrho }_2$	${\varrho }_3$	${\varpi }_1$	${\varpi }_2$	${\varpi }_3$
0.04	0.01	0.05	0.05	0	0	0.05	0	0

Table 2. Table 2 : Parameters of the model used for the numerical experiment 3 .

$T,\text{yr}$	$E,\$$	$\gamma$	$\kappa$	$\overline{v}$	${\eta }_d$	${\eta }_f$	$r_{d,0}$	$r_{f,0}$
0.25	100	0.3	0.5	0.1	0.007	0.012	0.1	0.1
$v_0$	${\lambda }_d$	${\lambda }_f$	${\varrho }_1$	${\varrho }_2$	${\varrho }_3$	${\varpi }_1$	${\varpi }_2$	${\varpi }_3$
0.04	0.01	0.05	0.074	0.014	2.10	1.0	0.5	0.5

Table 3. Table 3 : The ROC in the s 𝑠 s direction observed in the numerical experiment ( 1 ) by using PM.

$m_1$	$h_{\min }$	$c_s$	$V_1$	ROC	$V_2$	ROC
8	13.02	1834.59	8.25786	-	7.73678	-
16	5.84	1130.55	8.45893	-	7.92453	-
32	2.78	627.29	8.42466	2.55	7.89247	2.55
64	1.36	330.29	8.41957	2.75	7.88750	2.68
128	0.67	169.45	8.42030	2.79	7.88808	3.09
Mean of ROC				2.7		2.7

Table 4. Table 4 : The RE ϵ italic-ϵ \epsilon obtained in the experiment 3 as a function of the number of nodes in various directions, and the elapsed time.

$m_1$	$m_2$	$m_3$	$m_4$	$\Delta \tau$	$V_1$	$V_2$	${ϵ}_1$	${ϵ}_2$	El. time
8	6	6	6	0.01	4.120	4.048	3.04E-2	3.05E-2	3.21
10	8	8	8	0.005	3.746	3.681	6.30E-2	6.28E-2	19.11
12	10	10	10	0.0025	3.805	3.738	4.84E-2	4.83E-2	88.77
16	14	10	10	0.002	3.975	3.906	5.88E-3	5.82E-3	256.30
20	14	10	10	0.000625	4.006	3.936	1.75E-3	1.80E-3	647.73

Table 5. Table 5 : The RE ϵ italic-ϵ \epsilon obtained in the experiment 1 as a function of the number of nodes in various directions, and the elapsed time.

Method	$m_1$	$m_2$	$m_3$	$m_4$	$\text{Re}\left({\lambda }_{\max }\right)$	$V_1$	$V_2$	${ϵ}_1$	${ϵ}_2$	El. time
FDKM
	10	8	6	6	-626.43	20.552	19.136	1.44E-0	1.42E-0	0.22
	20	16	12	12	-4361.31	13.903	12.966	6.51E-1	6.43E-1	9.19
	20	16	14	14	-4361.36	13.904	12.966	6.51E-1	6.43E-1	21.17
	28	20	14	14	-9913.52	5.758	5.409	3.16E-1	3.14E-1	49.64
	34	24	20	20	-15670.10	9.811	9.170	1.65E-1	1.62E-1	320.78
	36	28	20	20	-17967.50	10.157	9.491	2.06E-1	2.03E-1	446.20
PM
	10	8	6	6	-296.29	8.300	7.785	1.42E-2	1.30E-2	0.29
	20	16	12	12	-4184.96	8.444	7.910	12.86E-3	2.80E-3	7.91
	20	16	14	14	-4200.14	8.443	7.910	2.72E-3	2.80E-3	11.06
	28	20	14	14	-10567.70	8.438	7.905	2.15E-3	2.19E-3	38.30
	34	24	20	20	-17827.90	8.437	7.903	2.09E-3	1.94E-3	275.44

Table 6. Table 7 : The RE ϵ italic-ϵ \epsilon obtained in the experiment 2 as a function of the number of nodes in various directions, and the elapsed time. < 0 absent 0 <0 means the price is negative.

Method	$m_1$	$m_2$	$m_3$	$m_4$	$\text{Re}\left({\lambda }_{\max }\right)$	$V_1$	$V_2$	${ϵ}_1$	${ϵ}_2$	El. time
FDKM
	10	8	6	6	-			-	-	-
	16	10	8	8	-			-	-	-
	20	16	12	12	-4495.91	13.951	12.875	1.13E-1	2.15E-1	14.34
	24	18	14	14	-7056.19	12.453	11.484	6.03E-3	8.39E-2	37.82
	28	20	16	16	-10211.90	5.806	5.319	5.36E-1	4.97E-1	115.69
PM
	10	8	6	6	-1926.40	12.498	10.601	2.39E-3	6.84E-4	0.40
	16	10	8	8	-8104.68	12.500	10.582	2.26E-3	1.14E-3	3.18
	20	16	12	12	-14529.30	12.548	10.612	1.58E-3	1.73E-3	20.33
	24	18	14	14	-22735.3	12.519	10.590	7.86E-4	3.81E-4	57.50
	28	20	16	16	-32725.50	12.533	10.598	3.78E-4	3.98E-4	171.01

Equations

$$\begin{array}{rl}d{r}_{t,d}=& {\lambda }_d({\theta }_d(t)-r_{t,d})dt+{\eta }_{d}d{W}_{t,d}^{\mathbb{Q}},\\ d{r}_{t,f}=& {\lambda }_f({\theta }_f(t)-r_{t,f})dt+{\eta }_{f}d{W}_{t,f}^{\mathbb{Z}},\end{array}$$

$$\frac{d{s}_t}{s_t}$$

$$d{v}_t$$

$$d{r}_{t,d}$$

$$d{r}_{t,f}$$

$$s_t\in [0,\infty ),r_{t,d}\in (-\infty ,\infty ),r_{t,f}\in (-\infty ,\infty ),v\in [0,\infty ),t\in [0,\infty ).$$

$$\displaystyle\langle d\mathbf{W}_{t}d\mathbf{W}^{\top}_{t}\rangle=\left(\begin{array}[]{cccc}1&\rho_{s,v}&\rho_{s,d}&\rho_{s,f}\\ \rho_{s,v}&1&\rho_{v,d}&\rho_{v,f}\\ \rho_{s,d}&\rho_{v,d}&1&\rho_{d,f}\\ \rho_{s,f}&\rho_{v,f}&\rho_{d,f}&1\\ \end{array}\right)dt,$$

$$\frac{\partial V}{\partial \tau }=$$

$$\mathcal{L}V=$$

$$+{\rho }_{s,v}\gamma sv\frac{{\partial }^{2}V}{\partial s\partial v}+{\rho }_{s,d}{\eta }_{d}s\sqrt{v}\frac{{\partial }^{2}V}{\partial s\partial r_d}+{\rho }_{s,f}{\eta }_{f}s\sqrt{v}\frac{{\partial }^{2}V}{\partial s\partial r_f}$$

$$+{\rho }_{v,d}\gamma {\eta }_d\sqrt{v}\frac{{\partial }^{2}V}{\partial v\partial r_d}+{\rho }_{v,f}\gamma {\eta }_f\sqrt{v}\frac{{\partial }^{2}V}{\partial v\partial r_f}+{\rho }_{d,f}{\eta }_d{\eta }_f\frac{{\partial }^{2}V}{\partial r_d\partial r_f}$$

$$+(r_d-r_f)s\frac{\partial V}{\partial s}+\kappa (\bar{v}-v)\frac{\partial V}{\partial v}+{\lambda }_d({\theta }_d(\tau )-r_d)\frac{\partial V}{\partial r_d}$$

$$+\left({\lambda }_f({\theta }_f(\tau )-r_f)-{\rho }_{s,f}{\eta }_f\sqrt{v}\right)\frac{\partial V}{\partial r_f}-r_{d}V,$$

$$V(0,s,v,r_d,r_f)={\left(s-E\right)}^{+},$$

$$V(0,s,v,r_d,r_f)={\left(E-s\right)}^{+},$$

$$[\tau ,s,v,r_d,r_f]$$

$$\equiv (0,T]\times \Omega ,$$

$$\varphi (\Vert \mathbf{x}-{\mathbf{x}}_i{\Vert }_2)=\exp \left[{\left(-\frac{\Vert \mathbf{x}-{\mathbf{x}}_i{\Vert }_2}{c}\right)}^2\right],i=1,2,\dots ,m,$$

$$\{x_i-h,x_i,x_i+{\omega }_{i+1}h\},{\omega }_{i+1}>0,h>0,$$

$$f^{\prime }(x_i)\simeq {\alpha }_{i-1}f(x_{i-1})+{\alpha }_{i}f(x_i)+{\alpha }_{i+1}f(x_{i+1})={\hat{f}}^{\prime }(x_i),2\le i\le m-1,$$

$${\alpha }_{i-1}$$

$${\alpha }_i$$

$${\alpha }_{i+1}$$

$$\begin{array}{rl}\{x_{i-2},x_{i-1},x_i,x_{i+1}\}& =\{x_i-w_{i-2}h,x_i-h,x_i,x_i+w_{i+1}h\},\\ & w_{i-2},w_{i+1},h>0.\end{array}$$

$$f^{\prime \prime }(x_i)={\beta }_{i-2}f(x_{i-2})+{\beta }_{i-1}f(x_{i-1})+{\beta }_{i}f(x_i)+{\beta }_{i+1}f(x_{i+1}),3\le i\le m-1.$$

$${\beta }_{i-2}$$

$$\displaystyle-h^{2}w_{i-2}\left(\left(w_{i+1}-3\right)w_{i+1}+1\right)\Big{\}}\Big{\{}c^{2}h^{2}\left(w_{i-2}-1\right)w_{i-2}\left(w_{i-2}+w_{i+1}\right)\Big{\}}^{-1},$$

$${\beta }_{i-1}$$

$${\beta }_i$$

$${\beta }_{i+1}$$

$$\displaystyle-h^{2}\left(w_{i-2}\left(w_{i-2}+3\right)+1\right)w_{i+1}\Big{\}}\Big{\{}c^{2}h^{2}w_{i+1}\left(w_{i+1}+1\right)\left(w_{i-2}+w_{i+1}\right)\Big{\}}^{-1},$$

$${\mu }_1$$

$$+h^2\left(w_{i-2}+1\right)w_{i+1}^2$$

$${\mu }_2$$

$$+h^2\left(w_{i+1}-1\right)w_{i-2}^2.$$

$$\begin{array}{rl}& c_s=2\max \{\mathbf{\Delta }\mathbf{s}\},c_v=3\max \{\mathbf{\Delta }\mathbf{v}\},\\ & c_{r_d}=3\max \{{\mathbf{\Delta }\mathbf{r}}_d\},c_{r_f}=3\max \{{\mathbf{\Delta }\mathbf{r}}_f\},\end{array}$$

$$M_{s}=\left(\alpha_{i,j}\right)_{{m}\times{m}}=\left\{\begin{array}[]{lll}\alpha_{i,j}\ \text{from}\ (\ref{1Dweight3unstructured1})&&\quad i=j,\\ \alpha_{i,j}\ \text{from}\ (\ref{1Dweight3unstructured2})&&\quad i-j=1,\\ \alpha_{i,j}\ \text{from}\ (\ref{1Dweight3unstructured3})&&\quad j-i=1,\\ 0&&\quad\text{otherwise}.\end{array}\right.$$

$${\alpha }_{1,1}={\alpha }_{m,m-1}=\frac{h}{c^2}-\frac{1}{h},{\alpha }_{1,2}={\alpha }_{m,m}=\frac{1}{h}.$$

$$M_{ss}=\left(\beta_{i,j}\right)_{{m}\times{m}}=\left\{\begin{array}[]{lll}\beta_{i,j}\ \text{from}\ (\ref{2Dweight3unstructured1})&&\quad i=j,\\ \beta_{i,j}\ \text{from}\ (\ref{2Dweight3unstructured2})&&\quad i-j=1,\\ \beta_{i,j}\ \text{from}\ (\ref{2Dweight3unstructured3})&&\quad j-i=1,\\ \beta_{i,j}\ \text{from}\ (\ref{2Dweight3unstructured4})&&\quad i-j=2,\\ 0&&\quad\text{otherwise}.\end{array}\right.$$

$${\beta }_{1,1}={\beta }_{m,m-1}=\frac{-4}{c^2},{\beta }_{1,2}={\beta }_{m,m}=\frac{2}{c^2}.$$

Full text

1 Introduction.

As per NASDAQ, the foreign exchange (FX) market is the most actively traded market in the world. More than $5 trillion are traded on average every day. By comparison, this volume exceeds global equities trading volumes by 25 times. Accordingly, the FX options market is the deepest, largest and most liquid market for options of any kind. Therefore, mathematical modeling of the FX options is an important area of modern mathematical finance. There exists a wide literature on this field, see, e.g. [32, 7] and references therein. From a risk management prospective, the FX models are useful, as they permit us to investigate the effects of severe market crashes on the FX rates. This is significant for long–dated (maturities of 20 years or more) FX derivatives embedded with popular early exercise contract features, see [15, 29] and references therein.

Having in mind that a majority of the FX options are exotic, nowadays it is a common approach for the FX models to treat the underlying domestic and foreign interest rates (IRs) to be stochastic as well as volatility of the FX rate itself, [23, 34, 43]. However, because of this complexity, pricing options under such a model requires using numerical methods, in more detail see [47, 26, 38] and references therein. While Monte Carlo methods are traditionally slow, the FD methods suffer from the curse of dimensionality and the FE approaches requires an extensive triangulization. Therefore, various attempts have been taken in the literature to propose a tractable yet sophisticated enough model to be able to capture an observed market dynamics of option prices.

For instance, the authors of [21] improved the Heston stochastic volatility (SV) model [22] in the FX setting under the postulate of constant domestic and foreign IRs. Although the assumption of constant IRs is definitely appealing due to its simplicity, empirical results have confirmed that such models do not reflect the market reality, specially in the case of a new generation of long–dated hybrid FX products. For these products, the fluctuations of both the exchange rate and the IRs are critical, so that the constant IRs assumption is clearly inappropriate for reliable valuation and hedging.

In [45] another improvement over the Heston SV model is presented. This model is then further extended in [39, 19] for currency derivatives by taking into account stochastic IRs and assuming all stochastic factors are correlated. A multi–factor SV model of the Heston type was introduced in [8]. Their model is coherent with respect to triangular relationships among currencies and allows for a simultaneous calibration of the volatility surfaces of the FX rates involved in a triangle, such as EUR/USD/JPY.

Despite these attempts, it was recognized that yet more sophisticated models could be helpful for modeling FX derivatives. In particular, there has been a great interest in modeling FX derivatives using four–factor jump–diffusion models, see [1] and the references therein. Typically, in these models, the spot FX rate and its variance follow a jump–extension of the Heston model [17], while the domestic and foreign IRs follow the one–factor Hull–White or Cox–Ingersoll–Ross (CIR) dynamics, [3, 25].

Therefore, based on this short survey, in this paper we consider a model proposed in [18]. However, while the authors of that paper tried to use some additional approximations to make it analytically tractable (via a closed–form solution for the characteristic function), in this paper we consider the whole model, also assuming non–zero correlation between all stochastic factors. Despite this model losses analytic tractability, it could be efficiently dealt with by using a novel RBF–FD numerical scheme of the second order.

Main contributions of this paper are as follows:

• •

To tackle efficiently solving a 4D time–dependent PDE, we propose an adaptive (non–uniform) discretizations in order to use as low the number of discretization points as possible while preserving a fixed accuracy. Our localized RBF–FD approach will arise in sparse matrices.

• •

By using non–equidistant stencils with three and four nodes for approximating the first and second derivative of the function, respectively, we prove that our RBF–FD approach acquires the second order of convergence in the internal points.

• •

A new strategy is proposed for the choice of the shape parameter. In particular, the shape parameters vary in each spatial variable (each dimension), and are also a function of the number of discretization points in the corresponding dimension.

• •

Approximation of the boundary conditions for the PDE is proposed to speed up computations.

• •

In case parameters of the model are already known, and calibration to, e.g., a term structure of swaptions is not required, a simple yet effective idea is imposed to approximate the time–independent coefficients by constant values, that subsequently allows the set of discretized ordinary differential equations (ODEs) to have a constant system matrix.

As compared, e.g., with a recent paper [36] (further MS) which deals with the numerical pricing of financial derivatives using RBF–FD method of high order and also with the non–uniform node layout, the differences are as follows:

• •

MS considers a 2D problem while here we attack a 4D problem.

• •

No closed–form formulas with for the weights are given in MS, while here we present closed–form expressions for the weights provided the second order approximations.

• •

The method of MS tries to take into account as many nodes in the neighborhood of a point as possible. Here we do not do this.

• •

In our method we derive the diagonal elements of the differentiation matrices in closed form. Hence, when the number of (uniform or non–uniform) nodes increases, there is no need building up the system (7) of the MS paper each time.

The remaining parts of this work are organized as follows. In Section 2 we describe the model and provide a PDE, so prices of FX vanilla options solve it under the appropriate boundary and initial conditions. In Section 3, special weights of the Gaussian RBF–FD method for the first and second derivatives are proposed. It is proved, that with these novel weights the scheme achieves the second order of convergence. Then, the numerical solution of the PDE is discussed in detail in Section 4. An analysis of the convergence and stability of this solution is also provided. Numerical discussions and reports are presented in Section 5. Through the computational aspects, we depict the accuracy of the new numerical procedure for pricing under SV and IRs, and investigate the influence of specific model parameters. Some asymptotic solutions are also discussed there. Section 6 concludes and outlines possible future works.

2 Model.

As mentioned, in this section we deal with a four–factor model where all the processes are represented by Itô’s diffusion with drift. An extended discussion on four–factor models and their usefulness in practice can be found, e.g., in [2].

We consider domestic and foreign IR processes, $r_{t,d}$ and $r_{t,f}$ which follow a Hull–White (short-rate) dynamics, [25] defined under their corresponding spot measures ($\mathbb{Q}$–domestic and $\mathbb{Z}$–foreign), respectively:

[TABLE]

where $W_{t,d}^{\mathbb{Q}}$ and $W_{t,f}^{\mathbb{Z}}$ are Brownian motions under $\mathbb{Q}$ and $\mathbb{Z}$, respectively. Parameters $\lambda_{d}$, $\lambda_{f}$ determine the speed of mean reversion to the mean-reversion levels $\theta_{d}(t)$, $\theta_{f}(t)$, and parameters $\eta_{d}$, $\eta_{f}$ are the volatility of volatility, or vol–of–vol.

Consider the spot FX rate, $s_{t}$, which is expressed in units of domestic currency, per unit of a foreign currency. We follow [18] where the dynamics for the interest rates in (1) is combined with the Heston model for $s_{t}$ (so the whole model is called FX–HHW) assuming non-zero correlations among all stochastic factors

[TABLE]

Here all stochastic processes are defined under the domestic risk–neutral measure, $\mathbb{Q}$, and $\gamma$ is the volatility–of–volatility parameter for the process $v_{t}$. Under the domestic–spot measure the drift of $r_{t,f}$ contains an additional term $-\eta_{f}\rho_{s,f}\sqrt{v_{t}}$, see [18]. The full correlation matrix associated with the dynamics (2) reads

[TABLE]

where $\mathbf{W}_{t}=[W_{t,s}^{\mathbb{Q}},W_{t,v}^{\mathbb{Q}},W_{t,d}^{\mathbb{Q}},W_{t,f}^{\mathbb{Q}}]^{*}$, and all elements $\rho_{i,j},\ i,j\in[s,d,f,v]$ of the correlation matrix are constant. The initial values of all processes in (2): $s_{0},r_{0,d}\equiv r_{d,0},r_{0,f}\equiv r_{f,0},v_{0}$, are parameters of the model. To guarantee a non–negativity of the instantaneous variance, a standard Feller condition is assumed to be satisfied, i.e. $2\kappa\theta/\gamma^{2}>1$, [31].

Following a standard no–arbitrage argument, the pricing PDE could be derived for the FX–HHW model and reads, [18, 26]:

[TABLE]

[TABLE]

where $V=V(\tau,s,v,r_{d},r_{f})$, $\tau=T-t$ is the backward time, and $\mathcal{L}$ is a linear differential operator of the second order.

The PDE in (8) should be solved subject to the initial and boundary conditions. The initial condition for the plain vanilla Call option reads:

[TABLE]

and for the plain vanilla Put option is:

[TABLE]

where $E$ is the option strike. The boundary conditions are discussed in more detail in Section 4.1.

3 Numerical solution of (8)

The multi–dimensional 4D PDE (8) includes non–affine terms, i.e., square roots and products. Therefore, solving it requires some numerical method (see e.g. [20, 46]), like FD, or meshfree RBF method. At the same time, the presence of four spatial variables for this time–dependent problem as well as the existence of six mixed derivative terms implies that any numerical method is supposed to be expensive. For instance, one can consider spatial discretization of (8) by using a FD approach (though this has not yet been studied in the literature) or a FE method, for example [55]. But for four–dimensional problems these methods already suffer from the high computational burden (e.g., due to curse of dimensionality).

Last years RBF schemes have gained a lot of interest in computational finance, [24, 40], because with RBF high resolution schemes could be constructed by using only few discretization nodes. This is especially helpful when solving various multi–dimensional problems, e.g., for models whose settings use several stochastic factors.

Meshfree methods (specifically, the localized ones) constructed based on an RBF approximation have been shown to perform better than standard FD methods for option pricing problems in one or more spatial dimensions, see e.g. [9] and the references therein. The main difficulty when using RBF collocation approach is the necessity to invert an ill–conditioned matrix arising due to a global RBF support. Authors in [53] discussed another approach using RBF for solving PDEs, that is by constructing locally supported operators approximating derivatives in the same manner as in the case of traditional FD scheme. Their approach is called RBF–FD method. This method is a local method resulting in a sparse linear system in contrast to the global RBF–schemes, which lead to ill–conditioned dense matrix systems, [37]. Thus, localized method can give not only high accuracies inherited from the use of RBFs but also sparse structures making them efficient for solving high-dimensional problems.

3.1 Localization.

We remind that (8) is defined on an unbounded domain $[\tau,s,v,r_{d},r_{f}]\in(0,T]\times[0,\infty)^{2}\times(-\infty,+\infty)^{2}$ subject to the initial condition $g(\cdot)$ as it is defined in (9)–(10). To solve (8) numerically, we truncate this unbounded domain to a finite domain as follows:

[TABLE]

where sub-indexes $\min,\max$ mark the minimum and maximum values chosen appropriately.

In [30] it is shown that the price of a European vanilla option (e.g., calculated by solving (8)) is estimated exponentially well by that of the associated to the barrier option (in $\log H$ where $H$ is the barrier), and, hence, solution of the problem defined on the truncated domain still exists. To be more precise, proposition 3.1 in [30] reveals that a continuously monitored barrier option price estimates that one of a European option arbitrarily well by extending the log–barrier. Since, a truncated domain is needed to calculate a barrier option price, in the present case of European options, truncating the domain is efficient once the truncation boundary is far enough from the points of interest.

We emphasize that choosing, e.g., the upper truncated boundaries $s_{\max}$, $v_{\max}$, ${r_{d}}_{\max}=-{r_{d}}_{\min}$, ${r_{f}}_{\max}=-{r_{f}}_{\min}$, sufficiently far away, reduces the error of moving the boundary conditions from the original boundary to the artificial boundary. However, on contrary larger computational domain needs a larger discretization width. Accordingly, this increases the error of the approximation of derivatives, or it demands for a large number of discretization points to get a required accuracy. A non–equidistant distribution of computational nodes could overcome this issue while making the discretization and computer coding more difficult. This will be discussed in more detail in Section 4.

3.2 Numerical analysis of the Gaussian RBF–FD weights.

While there exist various popular choices of the basis functions, see e.g., [11], in this paper we deal just with the Gaussian RBF. Investigation of other choices could be done in a similar way. The motivation of choosing the Gaussian RBF is partly because in the model (2) the interest rates have a marginal normal distribution. This, however, is less important for the localized version of the RBF method.

The Gaussian RBF is defined as, [5]:

[TABLE]

where $c$ is the shape parameter, $\bf{x}$ is the vector of coordinates representing a point in the (four–dimensional) space, $\|\cdot\|_{2}$ is the $L_{2}$ norm in this space, and ${\bf x}_{i},\ i\in[1,m]$ is a set of discretization (or so–called collocation) nodes.

According to the RBF–FD method, we need to approximate derivatives in (8) on a given stencil which could be constructed by using a subset of the discretization nodes. For doing that a certain weight is assigned to each node in the stencil. In turn, to find the weights of the RBF–FD formulas (for the first derivative) we consider a three point stencil in one dimension with $x$ being the corresponding coordinate:

[TABLE]

where $\omega_{i}$ is the corresponding $i$–th weight, and $h$ is a constant step. Obviously, (13) is just another way of representing a non–uniform grid. In case of polynomial functions used by the standard FD methods, approximations of the first and second order derivatives as well as mixed derivatives on this grid are well–known, [5, 27].

As applied to the RBF–FD method, the three–point non–uniform stencil $[x_{i}-h_{i},x_{i},x_{i}+h_{i+1}]$ can be first represented as in (13). Then the approximations of derivatives with the necessary order could be derived given the explicit form of the RBF to determine weights $\omega_{i}$ once $h\ll c$. Below we provide these expressions in the explicit form.

Theorem 1**.**

Consider a three–point approximation of the first derivative,

[TABLE]

where $m$ is the total number of nodes in this dimension, and $i$ runs across all internal nodes (so excluding the boundary nodes). For sufficiently smooth function $f$ and Gaussian RBF in (12) this approximation is of quadratic convergence rate if the weights are defined as

[TABLE]

The proof can be found in Appendix A.

The three–point approximation in (14) gives rise to a sparse matrix of first derivatives across given nodes. That is similar to the standard FD method, but with the added advantage that the RBF–FD method can naturally handle scattered node layouts.

Unfortunately, for the second derivatives a three–point stencil yields just a first order approximation. This, however, can be resolved if we involve more points into the stencil. The idea is to consider four adjacent points (for the interior nodes) to resolve the above–mentioned problem.

Theorem 2**.**

Consider the following set of nodes:

[TABLE]

Also consider the following approximation of the second derivative of some function $f(x)$

[TABLE]

This approximation for a sufficiently smooth function $f$ is of quadratic convergence rate when the weights are defined as

[TABLE]

where

[TABLE]

The proof of this theorem is given in Appendix B.

3.3 Shape parameters.

For a localized RBF–FD scheme the value of the shape parameters plays less important role than for the global RBF methods. However, its sharp selection may help in getting a better accuracy of the approximations.

Here, we propose to choose the shape parameters used for the calculation of weights in Theorems 1,2 in an adaptive way. This means that the shape parameters differ for each spatial variable (each dimension), and are also a function of the number of discretization points in the corresponding dimension. Our particular form of this dependence is as follows

[TABLE]

where $\bm{\Delta s}$, $\bm{\Delta v}$, $\bm{\Delta r}_{d}$ and $\bm{\Delta r}_{f}$ are the vectors of increments along $s$, $v$, $r_{d}$ and $r_{f}$, respectively. This way of choosing the shape parameter has an advantage over the classical ways where the shape parameter is either fixed or is chosen based on the condition number of an interpolation matrix. The values of $c_{s},c_{c},c_{r_{d}},c_{r_{f}}$ in (24) are chosen such to satisfy the condition of Theorems 1–2. At the same time we do not want them to be too large to eliminate losing precision. Also recall that when the shape parameters in (24) are high, e.g., once $c$ significantly differs from $h$, the RBF–FD approximations (15)–(17) & (20)–(23) tend to the FD formulas defined on non–uniform grids, [27].

4 Discretizations.

To solve (8) numerically we need to discretize it in both the time and space domains. For the temporal discretization we rely on the method of lines (MOL), [44]. The weights introduced for the Gaussian RBF–FD scheme in Section 3 can be used for the spatial discretization of (8). Our idea of discretizing adaptively, i.e., by focusing on the important areas of the solution, is not only to have the desirable accuracy at certain problematic areas of the solution, but also to reduce the number of the necessary discretization nodes, so subsequently handling a system of discretized equations of a moderate size.

Discretization of each differential operator in (8) yields a sparse differentiation matrix (DM). For the convection terms (the first derivatives) this matrix takes the form

[TABLE]

For the 4D problem such discretization is done in a systematic way using natural ordering, so to gain much speed. Each point of the 4D mesh corresponds to one row of the DM. In this way, all the weighting coefficients are gathered up into one sparse (banded) matrix.

The weights (15)–(17) can be used for the rows $2,3,\ldots,m-1$. If the boundary conditions are provided in the Dirichlet form, then there is no need for the first and last rows. If, however, we have a Neumann boundary condition (which is a reasonable choice for the interest rates), then for the first and last rows of the DM (25) the weights should be reconstructed as in this case the stencil will include ghost points outside of the chosen set of computational nodes. For instance, here with the RBF-FD methodology in use only two nodes are utilized instead of three. Doing in such a manner, we obtain

[TABLE]

This approach preserves a three–diagonal structure of the DM, but rigorously speaking breaks the second order of approximation at the boundary. Alternatively, one can replace the central approximation at the boundary with a one–side three–points approximation, because then the matrix still remains sparse (while not three–diagonal, however).

Using the result of Theorem 2 for the diffusion terms (the second derivatives), we obtain the following DM

[TABLE]

Again, in case of the Neumann boundary condition, this scheme for the first and last rows should be altered. In doing so, in this paper we again use an RBF–FD approximation

[TABLE]

Moreover, the discretization in (28) should also be changed for the second row of the DM as the four–point stencil now includes a ghost point. For instance, we can use three nodes instead of four to keep the diagonal structure of the DM, but loosing the second order of approximation, or use a one–sided approximation. The former approach based on the stencil $\{x_{1}-w_{1}h,x_{2},x_{3}+h\}$ gives rise to the following expressions

[TABLE]

The DMs produced based on (27)–(29) are also sparse (banded).

Spatial discretization of the cross derivative terms in (8) could be done by using the Kronecker product of the DMs, [33]. Therefore, in this way the weights for the corresponding approximations could be found by using (25) as a building block.

Gathering the weights coming from the convection, diffusion and source terms of the operator $\mathcal{L}$ in (8) all together results to the following system of ODEs

[TABLE]

with $A(\tau)$ being the entire matrix (sparse) containing rows and columns generated by all spatial dimensions. The weights used to construct this matrix allows the accuracy of the method to be similar to that of the global RBF methods. At the same time the sparse structure of $A(\tau)$ significantly reduces the computational cost of solving (30).

Also, it is worth mentioning that for the European options considered in this paper the payoff functions in (9), (10) have a discontinuity in the first derivative which prohibits high order convergence. In [36] these functions were smoothed using an established technique for Cartesian grids, namely by using a fourth order smoothing operator defined via a Fourier transform. However, here we do not use smoothing.

4.1 Boundary conditions.

For the European vanilla options, the usual boundary conditions in the $s$ space domain are

[TABLE]

Also, as the option price is usually not known at the boundary for the interest rates, it is a standard practice to set constant fluxes at the lower and upper boundaries, i.e., use a homogeneous Neumann boundary condition

[TABLE]

It is also well–known that if the Feller condition is satisfied, no boundary condition is required at $v=0$, the PDE itself with the value $v=0$ substituted into it should be used as the boundary condition, see e.g., [28] and references therein. However, if the Feller condition is not satisfied, the boundary condition must be set at $v=0$. That is usually either a reflection or a killing condition, [6].

At $v=v_{\max}$ there are various approaches to setting the boundary condition. One way is to proceed similar to (32) and set

[TABLE]

Alternatively, in [20] the authors propose to use

[TABLE]

which means that $V_{v}(\tau,s,v,r_{d},r_{f})=0$ at $v=v_{\max}$. In our experience, (33) provides more accurate results as compared with (34), for instance for the Heston model. That is because to be sufficiently accurate (34) requires the truncated boundary to be closer to infinity, while (33) relaxes this requirement.

By imposing the boundaries as in [49], we obtain the following system of linear homogenous (coupled) ODEs

[TABLE]

which should be solved subject to the initial condition given by a non–smooth payoff function in (9) or (10). Here $\bar{A}(\tau)$ is the coefficient matrix which includes the boundaries and is time–dependent. Note that it is singular as the boundary conditions are imposed directly into this matrix. Obviously (35) satisfies the Lipschitz condition, and thus a unique solution of (35) exists and extends to the whole working interval.

The following Lemma claims that the solution of (35) is conditionally uniformly stable.

Lemma 3**.**

Suppose the time–dependent matrix $\bar{A}(\tau)$ in (35) is of the size $N\times N$, where $N=m_{1}\times m_{2}\times m_{3}\times m_{4}$. Let us denote the largest and smallest point–wise eigenvalues of

[TABLE]

by $\lambda_{\max}(\tau)$ and $\lambda_{\min}(\tau)$. The solution of (35) is uniformly stable if there exists a finite constant $\delta$, such that the largest point–wise eigenvalue of (36) satisfies:

[TABLE]

for all $\tau$, $\tau_{0}$, such that $0\leq\tau_{0}\leq\tau$.

The proof reads a similar spirit of logic as in [42, page 133].

4.2 Krylov algorithm for the time–independent case.

Exponential time integration (ETI) schemes are schemes involving the matrix exponential and related matrix functions, [10]. An attractive feature of these integrators is a combination of excellent stability and accuracy properties, with the latter being usually better than in the standard explicit/implicit time integrators, see [13]. The interest in the ETI is due to the efficient programming of the Krylov subspace techniques to compute actions of matrix functions for large matrices.

As long as the coupled system (35) is time–independent, one may construct the solution by formally integrating (35), [41]:

[TABLE]

Since $\bar{A}$ is very large and sparse, a quick way of computing the solution is to rely on the Krylov subspace method. Further, without loss of generality, let $\tau=1$. Let us approximate $V(\tau)$ in (38) by constructing an orthonormal basis of the Krylov space

[TABLE]

and taking into account the orthogonalization of Gram–Schmidt (i.e., the algorithm of Arnoldi), [50]. If $V$ is an $N\times Y$ matrix with columns $v_{1},v_{2},\ldots,v_{Y}$ (the orthonormal basis vectors of $\mathcal{K}_{Y}$), it is possible to write

[TABLE]

Here $h_{Y+1,Y}=\|V_{Y+1}\|_{2}$, $Y$ is the dimension of the corresponding space, and the number of scalar products is proportional $Y^{2}$. Also $e^{*}_{Y}=(1,1,\ldots,1)_{1\times Y}$. Hence

[TABLE]

Here in fact, $H_{Y}$ is the restriction of $\bar{A}$ on $\mathcal{K}_{Y}$. Applying the Krylov approach makes it possible to project the main large scale problem onto the $\mathcal{K}_{Y}$ subspace and resolve it therein. As $v=\beta Ve_{1}$, $\beta=\|V(0)\|_{2}$, we have [52]:

[TABLE]

Considering $w_{Y}=\beta Ve^{H_{Y}}e_{1}$, we obtain

[TABLE]

where $VV^{*}w$ is the projection of $w=e^{\bar{A}}V(0)$ on $\mathcal{K}_{Y}$. In fact, $Y$ is much smaller than $N$ (the size of $\bar{A}$), and, therefore, the elapsed time required to calculate $e^{\bar{A}}V(0)$ is less.

Note that for time–stepping methods we need the largest eigenvalue of the system matrix to choose the best temporal step size, while in the described Krylov scheme this is not necessary. Also the described algorithm could be treated as computation of an action of a certain matrix function on the vector representing the initial condition. In more detail this algorithm is presented in Algorithm 1.

With Krylov approximations, considerable savings can be expected for large, moderately stiff systems of ODEs (35), which are generally solved by explicit time–stepping methods despite stability restrictions of the step size, or when implicit methods require prohibitively expensive Jacobians and linear algebra. It does not demand for any step size with no (strict) stability restrictions, see for example [12].

4.3 Integrating the time–dependent case.

Explicit time integration methods often require very small time steps if, e.g., the system of ODEs is stiff or a spatial mesh is locally refined at some points, [16]. However, when the system matrix is of a very large scale, this approach eliminates high computational burden inherent to implicit methods. In addition the explicit methods are easy to implement. Hence, we apply one of the most efficient methods of this type in the time–dependent case as it is described below in more detail.

Let us denote $V^{\iota}$ to be an approximate solution of the problem in contrast to the exact value $V(\tau^{\iota})$. Let us also consider $\varsigma+1$ equidistant temporal nodes with a temporal step size $\Delta\tau=T/\varsigma>0$, so $\tau^{\iota}=\iota\Delta\tau,\ 0\leq\iota\leq\varsigma$. To approximate the solution of the time–dependent problem (35) we use an explicit modified midpoint method by first applying the midpoint scheme [51]:

[TABLE]

Here $Z$ are intermediate approximations which march along in steps of $\Delta\tau$, and $V^{\iota}$ is the final approximation to $V(\tau^{\iota})$. This method is basically a centered difference or midpoint method except at the first and last points, and as such provides a second approximation in time. A simple yet efficient implementation of this scheme along its improved stability region can be found in [49, pages 119–123].

The motivation behind choosing this solver is to have a consistent second–order scheme in both space and time, which nowadays is a standard requirement in practice. Also (44) is a part of many standard mathematical software, e.g., in Wolfram Mathematica this can be done by the following call:

Method -> {"FixedStep", "StepSize" -> \Delta\tau,
      Method -> {"ExplicitModifiedMidpoint"}}

4.4 Choice of the time–dependent functions.

We define functions $\theta_{d}(\tau)$ and $\theta_{f}(\tau)$ in (8) as

[TABLE]

where $\varrho_{1}$, $\varrho_{2}$, $\varrho_{3}$ and $\varpi_{1}$, $\varpi_{2}$, $\varpi_{3}$ are six constants that can be found by calibration to the market data. However, using swaptions with different maturities for this purpose could result in a poor calibration. Therefore, an alternative is to replace (45) with a piecewise constant structure in time corresponding to the available maturities.

It is worth mentioning that our model has 16 parameters, therefore calibration of the model could be time consuming and even unstable. However, in the case when parameters of the model are already known, and calibration to, e.g., a term structure of swaptions is not required, a simple idea to improve performance of the model is to replace (46) assuming $\theta_{d}(\tau)$ and $\theta_{f}(\tau)$ are constants across the time. Our idea is to take the values of these constants to be the first term of the Taylor series expansion of (45) around $\tau=1$ and obtain

[TABLE]

Such an approach is not that practical, as again in such a form the model is not able to predict a term structure of, say swaptions. Therefore, in this paper this simplification is used only for testing the proposed method and investigation of its convergence and performance. Hopefully, in the general case of the time–dependent parameters $\theta_{d}(\tau)$ and $\theta_{f}(\tau)$ the characteristics of the method are assumed to be close to those found in this work.

5 Numerical experiments.

In this section we provide computational results obtained by using the proposed scheme (which further on is referred as PM). The code is implemented in Wolfram Mathematica 11.0 [14], and uses interpolation to find the required option value at any point of the space/time domain. In what follows we report solutions of (8) for the following instruments

A European vanilla Call option for parameters of the model given in Table 1 with $T=1$ year. A reference solution $V_{\text{ref}}(T,E,v_{0},0.024,0.024)$ $\simeq 8.420$ and $V_{\text{ref}}(T,E,v_{0},0.1,0.1)$ $\simeq 7.888$ has been obtained by using a very refined grid. 2. 2.

A European vanilla Put option for parameters of the model also given in Table 1, but with $T=2$ years. Here the reference values are: $V_{\text{ref}}(T,E,v_{0},0.024,$ $0.024)\simeq 12.528$ and $V_{\text{ref}}(T,E,v_{0},0.1,0.1)\simeq 10.594$. 3. 3.

A European vanilla Call option for parameters of the model given in Table 2. Here the reference prices are: $V_{\text{ref}}(T,E,v_{0},0.024,0.024)\simeq 3.999$ and $V_{\text{ref}}(T,E,$ $v_{0},0.1,0.1)\simeq 3.929$.

In all these numerical experiments the correlation structure of the model was defined as

[TABLE]

In each numerical experiment, in addition to the option price we also compute their Greeks. To remind, they are the derivatives of the option price on various model parameters, and are a very important tool for hedging an option portfolio risk associated with the market changes of these parameters, [54]. Calculation of Greeks is always a challenge for any numerical method, as a good approximation of the option price itself does not guarantee a good approximation for the Greeks.

Note that we did not compare our results with that in [56], since the latter includes only two and three factor models. As mentioned in Section 2, yet we are not aware of any numerical methods proposed in the literature for solving (8). Hence, for comparison, we also implemented an FD method. This method uses a three–point stencil and a central approximation of the first and second derivatives, and nine-point stencil for approximation of the mixed derivatives. It is known that the later approximation could be unstable and does not preserve positivity of the solution, see [27] and references therein, and we discuss this in more detail below in the paper when presenting the results. We use a uniform grid and the same time–stepping solver as in our RBF–FD scheme. The boundary conditions are also exactly same as that imposed for the PM. Overall, this scheme provides the second order of approximation in each spatial dimension, and the order of approximation in time depends on the time solver in use. We do not use any alternating direction implicit (ADI) approach, so all spatial nodes are combined into a single one-dimensional vector. This approach is chosen to exactly mimic what we do with the RBF–FD method. This implementation further on is referred as FDKM. We also set “$\mathtt{PrecisionGoal->5,AccuracyGoal->5}$” in our codes to speed up the processes as much as possible.

To investigate convergence properties of our scheme, we compute the standard relative error (RE) and the rate of convergence (ROC)

[TABLE]

where $V(m)$ is the solution of our problem obtained by using $m_{1}$ nodes. This expression assumes that the ROC is investigated separately for each spatial dimension.

All calculations were done under Windows 7 Ultimate, Intel(R) Core(TM) i5–2430M CPU 2.40GHz processor, HDD internal memory and 16.00 GB of RAM. The elapsed time is reported in seconds, and throughout the tables the notation $aE-b$ means $a\cdot 10^{-b}$.

The non–equidistant sets of nodes for $s,v,r_{d}$ and $r_{f}$ are defined as in [4, 20]

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

$m_{1},m_{2},m_{3},m_{4}\gg 1,1\leq i\leq m_{1},1\leq j\leq m_{2},1\leq k\leq m_{3},1\leq l\leq m_{4}$, and $x_{i},\ y_{j}$ are uniform points on $[0,1]$, see [48] for convergence properties of such non–uniform girds. The purpose of this kind of mesh generation is to have as many computational nodes as possible close to the most important (or problematic) areas. In particular, they include points close to the strike price (where the payoff function is not smooth); or at the lower boundary where the variance vanishes; or at the points corresponding to initial values of the domestic and foreign interest rates.

In our experiments we choose the computational boundaries to be set at $s_{\min}=v_{\min}=0,s_{\max}=14E,v_{\max}=10,r_{d,\min}=r_{f,\min}=-1,r_{d,\max}=r_{f,\max}=1$. When generating the nodes according to (49) we use $\xi_{s}=0.01,\xi_{v}=50,\xi_{rd}=\xi_{rf}=500$.

When the RBF–FD methods are in use, there is some freedom in placement of the computational nodes. Also in the areas close to the boundaries the nearest neighbor based stencils deform automatically, hence requiring no special treatment for computing the differentiation weights in those areas. However, from the computational efficiency prospective there exists a common intuition that one can use less discretization nodes in the $r_{d}$ and $r_{f}$ directions than in the $v$ direction, and similarly fewer discretization nodes in the $v$ directions than in the $s$ direction.

As suggested in [35], an alternative way to impose (31)–(33) in order to save computational time is as follows. One can consider the original PDE for the collocation nodes located at the boundaries, and discretize it by using the Gaussian RBF–FD approach. Then such discretized equations could be used themselves as the approximated boundary conditions to the PDE. In what follows we refer to this approach as the ABC conditions.

Below in our numerical experiments which deal with European vanilla Call options (those in 1 and 3), we impose only the Dirichlet boundary conditions. In contrast, for the European vanilla Put option in 2 for comparison reasons for all the nodes located at the boundaries we use the ABC conditions. As could be seen in Tables 6-6 and Figs. 1–2, this approach is computationally efficient from the performance point of view while preserving almost same accuracy of the solution.

5.1 Results.

Here, we discuss the results obtained in the numerical experiments 1–3.

In Table 6 the RE $\epsilon$ computed versus the reference value is presented for the FDKM and PM methods. Here $V_{1}=V(T,E,v_{0},0.024,$ $0.024),V_{2}=V(T,E,v_{0},r_{d0},r_{f0})$, and $\epsilon_{1},\epsilon_{2}$ are the corresponding REs. And in Table 3 the ROC of the PM method is given for several tests which differ by the number of nodes $m_{1}$ only in the $s$ dimension. The second order convergence (actually 2.7) can be observed.

The PM method does not include any special consideration to preserve positivity of the solution, as this was done, e.g., in [26]. Therefore, to verify it our results are also presented in Fig. 1 along with the graphs of Delta, Vega and Vanna all showing a stable behavior of the numerical solution.

For the European vanilla Put option 2, the results are brought forward in Table 7. As in this case the option maturity is longer, this clearly affects the largest eigenvalues of the system matrix. Note, that in the experiments 1–2 we are dealing with the constant system matrix in (35), and thus the Krylov scheme described in Algorithm 1 could be applied to calculate the solution. Numerical results obtained in this experiment are also presented in Fig. 2 to demonstrate the usefulness of PM and the way the artificial boundaries for our model are imposed.

Comparison of the results for test 2 with those obtained by using FDKM is furnished in Figure 3. For both methods the Krylov subspace method is used as the time–stepping solver. It can be seen that FDKM demonstrates some instability and even oscillations close to $v_{\max}$, while the RBF–FD method provides a stable solution in this area. A possible reason for these oscillations could be a nine–point approximation of the mixed derivative terms which does not preserve positivity of the option price.

[FIGURE:]

[FIGURE:]

The last experiment 3 is of importance since it gives rise to the time–dependent system matrix which makes the entire procedure of getting the numerical solution more difficult. The fast Krylov subspace method is no longer applicable here, and we rely on the explicit improved mid–point scheme to advance along time. Since, the FDKM is quite slow, here we use only PM and report the results in Table 4. The behavior of the method convergence looks similar to the previous experiments though for very accurate solutions much CPU time is consumed.

We also provided a test where the constant mean-reversion levels are obtained as in (46). This approximation as compared with the time–dependent levels is shown in Fig. 4. In the case of constant levels no time–stepping solver is required, and once again the Krylov subspace method can be utilized to provide a rapid convergence within a practically reasonable interval of time. The convergence of the method in this case is displayed in Table 6, and a quicker convergence could be observed as compared with that in Table 4.

Since the Krylov method is “step–size free”, i.e., there is no restriction on the step size, the solution could be obtained even in one time step. However, for marching along the time by using (44) in each case we used the time step given in Table 4. Once $\Delta\tau$ gets sufficiently small, a second order approximation can be reached as expected.

6 Final remarks.

In this paper, we consider pricing foreign exchange options by using a four–factor model including stochastic volatility and stochastic interest rates. Since this model is Markovian, the option price could be found by solving a four–dimensional PDE, which is of a parabolic type. Due to the complexity of the model this PDE should be solved numerically. The main idea of the paper is to propose a meshless method that could efficiently tackle this problem providing second order of approximation in both space and time.

The spatial discretizations are constructed by using Gaussian RBFs. This is because for the global RBF method they provide similarity to the marginal distributions of the interest rates in the considered model. However, we are concentrated on the local Gaussian RBF–FD method to provide a better performance. In particular, we explicitly derive the weights of the method which for non–equidistant nodes possess quadratic convergence.

The temporal discretization is done by using a method of lines. Hence, this construction results in a homogenous coupled set of ODEs with the system matrix being in general time–dependent. In this case an exponential time integrator is used to solve the system. If, however, the system is time–independent, the Krylov subspace scheme is applied to provide a quick convergence. A merit of the Krylov subspace scheme is its versatility. In addition, an adaptive procedure for selecting the involved shape parameters for each stencil is described.

Numerical experiments where the prices of the European Call and Put options are computed by solving the PDE, demonstrate good performance of the method while also prove the second order of approximations. Despite no special treatment of the initial condition (the option payoff) was provided roughly, the numerical solutions are smooth enough. We also compute option Greeks that demonstrate stable and qualitative correct behavior.

A natural extension of this work would be to adapt the method for pricing American–style FX options. Furthermore, in case where the system of ODEs (35) is stiff, the convergence of the Krylov method is expected to be slow. In this situation a (tensor–type) preconditioning for computing the action of this matrix on the payoff vector could be desirable. These will be investigated elsewhere.

Disclosure statement.

No potential conflict of interest is reported by the authors.

Appendix A Proof of Theorem 1

We use (14), take into account an explicit representation of function $f$ as the Gaussian RBF defined in (12), and apply it at the points $x_{i-1}=x_{i}-h,x_{i},x_{i+1}=x_{i}+\omega_{i+1}h$. With that one may obtain the weighting coefficients in (15)–(17) analytically in the limit $h\ll c$. With more rigor, if $\mu=h/c$, then the terms with the order $\mathcal{O}(\mu)$ should be also $\mathcal{O}(h^{2})$. The latter means that $c=\mathcal{O}(1/h)$. In what follows we assume that this condition is always satisfied by an appropriate choice of $c$.

To demonstrate the convergence order, we expand all weighting coefficients in (15) into series around $h=0$ up to the order two which yields

[TABLE]

Substitution of these expressions into (14) results in the following representation of the error provided by the approximation (14)

[TABLE]

This proves a second order approximation of (14) when the weights are defined as in (15)–(17).

Appendix B Proof of Theorem 2

Similar to the proof of Theorem 1, one may obtain the weighting coefficients (20)–(23) analytically assuming $c=\mathcal{O}(1/h)$. Substituting the Gaussian RBFs defined in (12) into (19) we obtain the following set of linear equations:

[TABLE]

where

[TABLE]

The solution to the system (54) is given by (20)–(23). Expanding the RHS into Taylor series around $h=0$ up to the second order yields:

[TABLE]

Substituting these expressions into (19) we obtain the following representation of the approximation error

[TABLE]

This proves that the proposed weights provide a second order of approximation for the second derivative $f^{\prime\prime}(x_{i})$.