Numerical Solution of Heston-Hull-White Three-Dimensional PDE with a High Order FD Scheme

: A new numerical method for tackling the three-dimensional Heston–Hull–White partial differential equation (PDE) is proposed. This PDE has an application in pricing options when not only the asset price and the volatility but also the risk-free rate of interest are coming from stochastic nature. To solve this time-dependent three-dimensional PDE as efﬁciently as possible, high order adaptive ﬁnite difference (FD) methods are applied for the application of method of lines. It is derived that the new estimates have fourth order of convergence on non-uniform grids. In addition, it is proved that the overall procedure is conditionally time-stable. The results are upheld via several numerical tests.


Introduction
To model different types of derivatives in finance, a common approach is to investigate the connections of these factors to each other, formulated as a stochastic differential equation (SDEs). The factors could be the underlying asset, the volatility, domestic and foreign interest rates, etc., [1,2]. As such, the important action of pricing option under different payoffs can be modeled and simulated via the SDEs or their corresponding partial differential equation (PDE) formulation.
However, a frequently occurring issue is that whatever the model becomes complicated and more realistic, the procedure of having and representing its exact solution becomes harder, see, e.g., [3][4][5].
To discuss more and from the beginning, the classical model of Black-Scholes in pricing contracts does not cover and illustrate all the aspects of an option in a complete market, such as market risks, stochastic volatility (SV), and asymmetries seen in data of market, [6]. Some remedies to this well-known model are via non-lognormal hypothesis for a SDE, that indicates some modifications of the volatility and the underlying asset. We recall that Heston in [7] extended and improved the behavior of the Black-Scholes model by involving more risky factor into the model, i.e., by considering the volatility to be stochastic as well. Further discussions can be found at [6,8].
On the other hand, as long as the foreign exchange (FX) products are involved and a trader encounters a situation in which the interest rate is not anymore constant during the lifetime of an option, then investigating and proposing an improved model, having stochastic rates of interest, such as the power-reverse dual-currency and the Heston-Cox-Ingersoll-Ross (HCIR) problems (refer to [9] and the references therein for more background).
In pricing under call options, the (terminal/)initial condition is given by [11,12]: where the strike price is E. In a similar way, for a put option, it is given as follows: As discussed in [13,14], the fair pricing procedure should be carried out by computational schemes since the corresponding high-dimensional PDEs, constructed for such options, do not admit any analytical or semi-analytical solutions, see [15,16] for further background.

Novelties and Motivation
The contribution of this article reads in proposing a solution method via an un-equally spaced grid having a focus on the hot area in option pricing under the HHW PDE problem. Studying and coding multi dimensional problems with discretization methods while the grid of points are non-uniform is a challenging and intensive task, but could clearly increase the accuracy of the approximate solution by applying fewer numbers of grid nodes in contrasts to the uniform discretization. This reduces the size of the discretized problem and is useful in practice.
To this aim, (1) is tackled by employing high order fourth-order finite difference approximations. We apply fourth order discretizations on a stencil having five and six non-equidistant nodes. Derivation and construction of fourth-order compact FD method for HHW PDE is new and useful in practice.
In fact, the method-of-lines technique is considered to build a set of ODEs with time-varying system matrix. All the side conditions are imposed therein as well. Thence, a method to march along time for the set of ODEs is provided in Section 3 and it is analytically illustrated that the presented numerical procedure is conditionally time-stable when b is not changing by time.
Recalling that here adaptive FD formulas are constructed to hit some features simultaneously, viz., to be effective, results in sparse operators and being able to handling non-uniform grids.
Motivated by recent works in this field (see e.g., [17]), we aim at proposing higher order schemes for the HHW equation on non-uniform meshes so as to increase the accuracy of obtained option prices without increasing the computational load so much. The novelties and contributions of our work are given below:

•
We propose fourth-order adaptive discretizations for the spatial variables.

•
The beauty of our scheme is the use of non-uniform grid of nodes with an adaptation on the hotzone.

•
We provide a new stability bound for the resulting fully discretized set of equations when pricing under HHW PDE using high order discretization methods along the spatial as well as the temporal variables.

Grid Generation
The option pricing problem (1) is considered in the unbounded area wherein For tackling the financial model numerically, one can take into account the following domain [18]: wherein s max , v max , r max are three positive real constants and assumed to be large enough. Since the PDE model is coercive (sometimes called degenerate) at v = 0, its payoff is non-smooth at s = E, and the working domain has large width, thus it is requisite to use non-uniform meshes, at which the location of the nodes are not equally-spaced. This helps in producing results of higher accuracy with adapting to the hotzone of the problem.
Let {s i } m i=1 be a set of non-uniform nodes along s as follows [13,19]: where m > 1 and ξ min = ξ 1 < ξ 2 < · · · < ξ m = ξ max are m equi-distant points with the following characteristics: , wherein s min = 0. Here d 1 > 0 controls the density of the nodes around s = E. We also have: Throughout this work, we used the same value for The nodes along v, i.e., {v j } n j=1 are defined by: where d 2 > 0 gives the concentration around v = 0. In this work, we used d 2 = v max 500 , where v max = 10. In addition, ς j are equally spaced points given by: for any 1 ≤ j ≤ n. The non-uniform nodes along r are defined as follows: whereas d 3 = r max 500 is a positive parameter and r max = 1. We also have ζ k = (k − 1)∆ζ, . Note that denser mesh points in the important area could circumvent the problems happening in solving (1), like non-smoothness of payoffs (2) and (3) at s = E, and the degeneracy at v = r = 0. We state that a detailed study into possibly better choices for the involved parameters in mesh generating may be interesting, but this is beyond the scope of the current research. Furthermore, the non-smoothness arising in the payoff would ruin the convergence rate of most derivative approximation particularly on uniform meshes and due to this, the application of non-uniform nodes is indispensable for efficient numerical solution of (1).

Manuscript Organization
The remaining parts of this work are organized as follows. In Section 2, the weights of the FD scheme over non-uniform grids (here we also call adaptive grids with special emphasis on the hot zone) are derived to attain the higher rate of convergence four. Section 3 is devoted to the application of a sixth order Runge-Kutta time stepping method to advance along time when semi-discretize the HHW PDE. We prove that the new procedure is time-stable conditionally based on the largest eigenvalue of the system matrix. Section 4 shows that numerical performances are more useful than the earlier schemes with quicker convergence behavior. Finally, some conclusions are drawn in Section 5.

Calculating the Weights of the High Order FD Scheme
In this section, by applying a methodology as in ( [20], Chapters 3-4) or [21], but with more Taylor expansion terms, we can construct fourth-order FD approximations on (non-uniform) grids.
Five points are required in estimating the first derivative as well as six points in approximating the second derivative in order to obtain a consistent fourth-order scheme throughout the discretized mesh of points.
Without losing the generality, let us construct the weights in the one dimensional case. Then, the concept of tensors using Kronecker product may be applied easily to transfer the weights to the appropriate dimensions. To this objective, consider a sufficiently smooth function g(s) and a grid as follows: {s 1 , s 2 , · · · , s m−1 , s m }.
Consider the following five adjacent nodes: and calculate the interpolation polynomial p(z) going via the nodes and then its first derivative p (z). At this moment, by employing a computer algebra system to do some symbolic computations and setting z = s i , we attain the FD estimate for the first derivative as follows: where the maximum local grid spacing is h and we have using Γ l,q = s l − s q and Recalling that the above procedure should be similarly done for the nodes {s 1 , s 2 , s m−1 , s m }, viz, to find the weighting coefficients with fourth order of convergence for such nodes, we should consider the five adjacent points and then calculate the interpolating polynomial at that specific point. In this way, the sided FD formulas are constructed and used.
Similarly, FD estimates for the second derivative terms can be obtained applying a similar methodology as above. To this objective, we consider a set of points as follows: and compute the second-derivative interpolating polynomial p (z) based on z. Now by taking into account z = s i in Mathematica [22], one obtains that: where , , Here, we have . Summarizing the following theorem has been established.
Theorem 1. As long as the function g is sufficiently smooth, the first and second derivative of the this function can be approximated by five and six adjacent points respectively on non-uniform meshes, via the formulas (13) and (17).

Proof. The proof can be investigated by Taylor expansions as in the derivation in this section.
It is hence omitted.
The procedure for obtaining the weights for the points {s 1 , s 2 , s 3 , s m−1 , s m } to keep the fourth convergence order should be investigated by the six adjacent points as described above but for that specific node.
It is noted that the formulations derived in (13) and (17) can be used for both uniform and nonuniform distribution of the discretization nodes, and can be simplified to more simpler formulations if the nodes are equidistant.

Application to Option Pricing under 3D HHW PDE
Considering the non-uniform nodes discussed in Section 1 along with the high order FD formulations calculated in Section 2, one is able to derive the differentiation matrices corresponding to the first and second derivatives of the function. These derivative matrices contains the weights of the fourth order approximations and are sparse in general since they are banded matrices whose zero elements are much more than their non-zero elements. These feature would help us in solving the financial model (1) as would be observed later.
For multi-dimensional derivatives, a matrix is constructed such that this is done on the flattened data, and subsequently the Kronecker product of the matrices for the derivatives (in one-dimension) are being considered.
One way for imposing the impact of (13)- (17) is with matrices including the weights of (13)-(17), i.e., the non-equidistant second-order FD weights, as their elements. A matrix which shows an estimation to the differential operator is called as a matrix of differentiation [20]. Forming and implementing the proposed scheme based on these matrices are invaluable aids for analysis.
Taking all the weights into consideration, the PDE (1) can be semi discretized to obtain: at which U(t) = (u 1,1,1 (t), u 1,1,2 (t), . . . , u m,n,o−1 (t), u m,n,o (t) Here the boundaries along s are defined as follows [13]: For v = v max , the following Dirichlet condition is prescribed: Remarking that the nodes which are located on the boundary v = 0 are considered as interior nodes and we take a fact into consideration that they must read the PDE model. That is to say, we incorporate the semi-discretized equations at this boundary.
By incorporating the above mentioned conditions, we obtain the following system of semi-discretized ODEs as follows: whereĀ(t) is the coefficient matrix including the boundaries.

Integrator
For discretizing in temporal variable t, many schemes are existing, for example refer to [23]. Explicit methods are basically straightforward to implement, but suffer from stability problems. Implicit schemes are unconditionally stable, but only exhibit low convergence or very time-consuming because of solving nonlinear system of algebraic equations per step.
Consider u ι to be the computational solution for the exact solution U(t ι ) and choose k + 1 temporal nodes with the step size ∆t = T k . At the moment, we use the δ-stage Runge-Kutta scheme [23] at t ι+1 = t ι + ∆t, (0 ≤ ι ≤ k) by: wherein f is defined based on the right hand side of (26). It is generally assumed that the row-sum conditions hold: Now we consider a sixth-order explicit Runge-Kutta scheme (RK6) below [24]: with b = (9/180, 0, 64/180, 0, 49/180, 49/180, 9/180), C = (1, 1/2, 2/3, (7 − p 1 )/14, (7 + p 1 )/14, 1), and p 1 = 21 1/2 . Notice that a consequence of explicitness is c 1 = 0 in (28), so that the function is sampled at the beginning of the current integration step. Here, the sixth-order time-stepping solver consists of seven stages and reaches sixth order of convergence. The sixth order says that the error of local truncation is on the order of O(∆t 7 ), while the total accumulated error is on the order of O(∆t 6 ).
In the sequel, we study that under what criteria the numerical discretized solution does not blow up. The following theorem is one of the contributions of this work. This is given for the time-independent case, i.e., whenĀ(t) =Ā. Proof. To find a stability conditions, we proceed as follows. Incorporating the time-stepping solver (27) on the system of ODEs (26) yields: Thus, the numerical stability is reduced to: which is due to (30) for any ω i as the eigenvalue ofĀ. Now by considering: ω max (Ā) , the stability condition can now be represented as follows: Noting that the negative semi-definiteness ofĀ makes all its eigenvalues to have negative real parts. Thus, the proposed scheme has numerical stability if the temporal step size ∆t satisfy (32). Noting that this can be computed in the language Mathematica [22] via the command The proof is ended.

Experiments
In this section, some tests were given for our proposed method showed via Adaptive Finite Difference Method (AFDM) to price at the money call options, when T = 1 year and E = 100$. A comparison was done by the standard uniform FD scheme [4], which by second order FD approximations and the Euler's scheme as a temporal solver shown by FDM. We also compare with the method provided in [13] shown by Haentjens-In't Method (HIM).
Mathematica 11.0 is used for the simulations [25]. Time is also reported in second while we employ the following stopping condition: wherein u ref and u approx are the exact and numerical results.
To increase the computational efficiency for very large scale semi-discrete systems that we are dealing with, here we set AccuracyGoal → 5, PrecisionGoal → 5.
Here, we consider more number of discretization nodes along s rather than v and r, since its working interval is larger than the others and the non-smoothness of the initial condition occurs along this spatial variable.
The non-constant b is defined as follows: where c 1 , c 2 , c 3 are constants, and τ = T − t. The following two test cases are considered: The results are brought forward in Tables 1 and 2 showing the stable and efficient valuations of options under HHW PDE via the new high-order procedure. Furthermore, to reveal the positivity and stability of the numerical results, in Experiment 2, and by considering m = 30, n = 18 and o = 18 discretization nodes, the results based on AFDM are plotted in Figures 1 and 2.

Ending Comments
In financial engineering, it is famous that the Black-Scholes PDE could not be useful in real application due to several restrictions. Several ideas to observe the market's reality are models based upon the stochastic volatility and interest rate models. The resulted PDE problem in this way is hard to be solved theoretically due to higher involved dimensions and so numerical methods are required.
In this paper, we have proposed a new discretized numerical method based on adaptive FD methodology on non-uniform grids in order to tackle an important problem in computational finance known as HHW PDE (1). It was proved that the new procedure has conditional stability and shown to be efficient in practice.
Further discussions can be investigated to extend the results of this work for other types of options defined on HHW model such as digital (binary) options, at which the initial condition is not only non-smooth at the strike but also discontinues.
Funding: This research received no external funding.