An Exploration of a Balanced Up-downwind Scheme for Solving Heston Volatility Model Equations on Variable Grids

This paper studies an effective finite difference scheme for solving two-dimensional 1 Heston stochastic volatility option pricing model problems. A dynamically balanced 2 up-downwind strategy for approximating the cross-derivative is implemented and analyzed. 3 Semi-discretized and spatially nonuniform platforms are utilized. The numerical method 4 comprised is simple, straightforward with reliable first order overall approximations. The 5 spectral norm is used throughout the investigation and numerical stability is proven. Simulation 6 experiments are given to illustrate our results. 7


Introduction
Demands of highly effective, efficient and reliable numerical methods have been increasingly high for solving option trading modeling equations involving cross-derivative terms.However, desirable computational procedures are in general difficult to obtain due to challenges from the participation of cross-derivatives [6,15].This motivates our study.In this investigation, targeting at European options that can only be exercised on dates of maturity, we propose and analyze a new and dynamically balanced up-downwind finite finite difference method in the pursuit.
In the early 1970s, Black, Scholes and Merton introduced the popular Black-Scholes-Merton (BSM) model [2,5].Under the consideration, stock prices are assumed to follow geometric Brownian motion, while the volatility of the stock prices is fixed and no sudden jumps occur.
However, classic BSM models often cannot fit ideally into market data observed nowadays [5].
This may be due to the fact that, in modern financial markets, not only stock prices are subject to risk, but also the estimate of the riskiness is typically subject to significant uncertainty.To incorporate additional sources of randomness into an option pricing model, Heston proposed a different approach by introducing the consideration of stochastic volatility [9].
There have been numerous recent publications on the numerical solution of Heston modeling equations.For instance, certain up-downwind first order algorithms are proposed and studied by Ma and Forsyth [15].Stability analysis are also carried out via standard von Neumann analysis for Cauchy problems or problems with periodic boundary conditions [4,12].Difficulties for more general stability analysis are primely due to the use of cross-derivative and boundary data [13].
Consequently, there has been no rigorous mathematical proof of the numerical stability for any second order scheme.
But cross-derivatives are essential to partial differential equations modeling a Heston Process.
Further, Heston modeling formulations also require more realistic Dirichlet, Neumann, or mixed boundary conditions [1,9].These have motivated our approaches.In this paper, we are particularly interested in computations based on a Heston put option model [4,5,8,10,12,22].
Similar investigations can also be carried out for a call option.
In particular, we consider the following two-dimensional Heston volatility model interpreting the behavior of the asset value S and its volatility y at time t ≥ 0,

dS(t) S(t)
= µdt + y(t)dW 1 (t), where µ is the expected return of the asset, κ is the rate of reversion to the mean level of the volatility, η is the mean level of the volatility, σ > 0 is the volatility parameter, and cov(u, v) is the covariance between u and v [9,21].The two Wiener processes W 1 (t) and W 2 (t) describe the random noise in asset and volatility, respectively.They are assumed to be correlated with a Let v(S, y, t), t ≥ 0, denote the value of a European put option that is a function of asset price, volatility and time.An application of Itô's Lemma and the non-arbitrage principle with a construction of risk-less portfolio leads to [5,7,9,12,19], Let v(S, y, T) = max {K − S, 0} , S, y ≥ 0, be the terminal condition to use, where T is the payoff time and K is the strike price.We adopt the following mixed boundary conditions for S, y > 0 and T > t ≥ 0 [4]: Set τ = T − t.Equation (1.1) can be rewritten as together with constraints [4,5,22], ) We may extend the temporal domain for (1.2)-(1.7)by allowing T = ∞.Further, for the sake of computations, we consider a truncated spatial domain Ω = {(x, y) : −X < x < X; 0 < y < Y} for sufficiently large X and Y in the rest of our investigations.
In the next section, a nonuniform spacial mesh will be introduced.Based on it, a semi-discretized system will be derived for solving (1.2)-(1.7).Dynamically balanced up-downwind difference approximations will be presented.A general linear stability analysis will be implemented in Section 3. Computational experiments will be carried out in Section 4.
Computationally evaluated rates of convergence of the scheme will also be provided.Finally, conclusions and future research intentions will be given in Section 5. in the -direction, respectively, where ∈ {x, y} [13,20].Similarly, for appropriate indexes, we define

Balanced up-downwind semi-discretized scheme
We now approximate the diffusion terms in (1.2) by using the above, and derivatives in (1.6) and (1.7) via the following, We approximate the advection terms in (1.2) through three different channels depending upon relations between η and r.

.6)
Define We now approximate the cross-derivative in (1.For the smoothness of nonuniform grids [20], we require that We propose that Substitute all spacial derivative approximations into (1.2) and let w denote the approximate solution to u.We acquire the following linear system, where w, f ∈ R MN and A ∈ R MN×MN is block tridiagonal in the form of ,  where It is observed that in the event if ρ = −1, we have the following due to (2.7): which indicate that uniform spacial grids must be employed.Thus, (2.9) reduces to Nontrivial entries of A s are readily to obtain based on above discussions.We need the following restrictions on mesh steps in the case [20]: Apparently, when ρ = 1, the above implies that a uniform spacial mesh with h = σk must be used.
Different from (2.8), we consider a new dynamically balanced cross-derivative approximation, In this circumstance, we obtain the following new system, where w, f (τ) ∈ R MN and Ã ∈ R MN×MN is block tridiagonal, that is, .
Nontrivial entries of Pm , Dm and Qm within their respective ranges of m are given by The semi-discretized method (2.12) reduces to a uniform scheme when ρ = 1, that is, Nontrivial elements of Ã can be determined from simplifications of the above formulae.

Lemma 3. [17]
The matrix exponential, e tA , tends to a zero matrix as t → +∞ if and only if all the eigenvalues of A have strictly negative real parts.
Proof.We will only need to show the case of ρ ∈ (0, 1], η > 2r for (2.9), since extensions of our results for other cases are technically imminent.Thus, we only need to show that each of the MN Geršhgorin discs of A lies on the left side of the complex plane.In fact, there are five types of the Geršhgorin discs to consider: 1. discs centered at an internal mesh point; We provide detailed proofs for the first three types of discs.Similar arguments can be applied to the rest cases.
CASE 1: In this situation, we first consider the situation in which η < y n ≤ Y. Let z ∈ S i be any complex number, where S i is a Geršhgorin disc centered at an internal point of the spacial grids.Thus,

.2)
Let α be the real part of z.Since we are concerned only about the upper bound of the real part of the eigenvalues, we may replace z by α via a triangle inequality, and remove absolute value sign on the left hand side of (2.2).As a consequence, (2.2) renders to Recall (2.10) and that ρ > 0. We have The above indicates that Furthermore, since y > η > 2r, we conclude that Therefore, the term inside each pair of absolute signs in (2.3) must be positive.We may remove all absolute signs in (2.3), and, subsequently, yields which is what we expect.Generalizations of the discussion for cases involving y ≤ η are straightforward.Therefore all eigenvalues contained in S i must lie on the left half of the complex plane.Peer-reviewed version available at Algorithms 2019, 12, 30; doi:10.3390/a12020030CASE 2: Without loss of the generality, we consider the case of x = x 1 and η < y < Y. Thus, for any complex number z ∈ S i , where S i is a Geršhgorin disc satisfying

Preprints
Similar to the previous case, we take α, the real part of z.Thus, The above apparently implies that such an S i must lie strictly on the left half of the complex plane, and the origin cannot be on its boundary.This ensures our expectation.
CASE 3: In the circumstance, Geršhgorin discs, S i , concerned are centered at boundary points where a Neumann condition is imposed.Hence, for any z ∈ S i we have The above indicates that α, the real part of z, must satisfy Recall Lemma 3.2.Since the origin cannot lie on the boundary of every Geršhgorin disc, combining results from the three cases, we conclude immediately that all eigenvalues of A must be strictly on the left half complex plane.Thus, we must have lim The above completes our proof.
key parameter value used strike price K = 100 interest rate r = 0.05 mean reversion speed κ = 2 long-run mean of volatility η = 0.1 market, we proceed with ρ = −1 and T = 5.For demonstrating the numerical solution and its rate of convergence estimates, we first consider uniform spacial grids.To this end, we may denote that h m = h, k n = k = σh, m = 1, 2, . . ., M; n = 1, 2, . . ., N.
Results over nonuniform grids will be presented later on.
Some key parameters used are shown in Table 1.Further, a Crank-Nicolson type temporal integrator will be utilized for advancing our semi-discretized system (2.9), (2.12), with ∆τ as the temporal step [18].It has been known that λ = ∆τ/c 2 , where c = min {h, k} , play an effective role of the Courant number [14,16].We experiment with different values of λ varying from 0.5 to 1.
Our semi-discretized scheme is expected to be up to the first order in convergence in space.To numerically examine this through experiments, we employ a generalized Milne's device [13,20].Then, for a selected terminal time T, we denote the numerical solution at point Let h = 0.01 and σ = 1.For simplicity of notations, we use the same letter v for the approximate solution to (1.1).We show the solution v for ρ = −0.5 and ρ = −1 in Figure 3 and Figure 5 , respectively.To see more precisely solution profiles, we show corresponding contour maps next to the surfaces.It can be observed that the European put option price is a decreasing function of the stock price S.This coincides well with the financial theory that a put option price should have a negative correlation with the underline stock price [1,11].To examine further the delicate relationship between a put option price and its volatility, we plot an average numerical solution v(y, t) taken across different stock prices with respect to the volatility in Figure 7.The simulated computational result is exactly what we would expect, since a put option price should be positively correlated with the volatility [11,16].

2 .
discs centered on one of the Dirichlet boundaries; 3. discs centered on the Neumann boundary; 4. discs centered at one of the intersection mesh points of two Dirichlet boundaries; 5. discs centered at one of the intersection mesh points of one Dirichlet boundary and the Neumann boundary.

Table 1 .
Key parameter values for numerical simulations