{'fraction_non_alphanumeric': 0.0418533292420988, 'fraction_numerical': 0.014482970236268794, 'mean_word_length': 4.137295081967213, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
High-order compact finite difference schemes for option pricing in stochastic volatility models on non-uniform grids

We derive high-order compact finite difference schemes for option pricing in stochastic volatility models on non-uniform grids. The schemes are fourth-order accurate in space and second-order accurate in time for vanishing correlation. In our numerical study we obtain high-order numerical convergence also for non-zero correlation and non-smooth payoffs which are typical in option pricing. In all numerical experiments a comparative standard second-order discretisation is significantly outperformed. We conduct a numerical stability study which indicates unconditional stability of the scheme.


Introduction
Efficient pricing of financial derivatives, in particular options, is one of the major topics in financial mathematics. To be able to explain important effects which are present in real financial markets, e.g. the volatility smile (or skew) in option prices, so-called stochastic volatility models have been introduced over the last two decades. In contrast to the seminal paper of Black and Scholes [BS73] the underlying asset's volatility is not assumed to be constant, but is itself modelled by a stochastic diffusion process. These stochastic volatility models are typically based on a two-dimensional stochastic diffusion process with two Brownian motions with correlation ρ, i.e. dW (1) (t)dW (2) (t) = ρ dt. On a given filtered probability space for the stock price S = S(t) and the stochastic volatility σ = σ(t) one considers dS(t) =μS(t) dt + σ(t)S(t) dW (1) (t), dσ(t) = a(σ(t)) dt + b(σ(t)) dW (2) (t), whereμ is the drift of the stock, a(σ) and b(σ) are the drift and the diffusion coefficient of the stochastic volatility.
In (2)-(4), κ * , v, and θ * denote the mean reversion speed, the volatility of volatility, and the long-run mean of σ, respectively. For some models and under additional restrictions, closed form solutions to (1) can be obtained by Fourier methods (see, e.g. [Hes93,Dür09]). Another approach is to derive approximate analytic expressions, see, e.g. [BGM10] and the literature cited therein. In general, however, -even in the Heston model when the parameters are non constant-equation (1) has to be solved numerically. Moreover, many (so-called American) options feature an additional early exercise right. Then one has to solve a free boundary problem which consists of (1) and an early exercise constraint for the option price. Also for this problem one typically has to resort to numerical approximations.
In the mathematical literature, there are a number of papers considering numerical methods for option pricing in stochastic volatility models, i.e. for two spatial dimensions. Finite difference approaches that are used are often standard, low order methods (second order in space). Other approaches include finite element-finite volume [ZFV98], multigrid [CP99], sparse wavelet [HMS05], or spectral methods [ZK10].
Let us review some of the related finite difference literature. Different efficient methods for solving the American option pricing problem for the Heston model are compared in [IT08]. The article focusses on the treatment of the early exercise free boundary and uses a second order finite difference discretization. In [IHF10] different, low order ADI (alternating direction implicit) schemes are adapted to the Heston model to include the mixed spatial derivative term. While most of [TGB08] focusses on high-order compact scheme for the standard (one-dimensional) case, in a short remark [TGB08, Section 5] also the stochastic volatility (two-dimensional) case is considered. However, the final scheme is of second order only due to the low order approximation of the cross diffusion term.
High-order finite difference schemes (fourth order in space) were proposed for option pricing with deterministic (or constant) volatility, i.e. in one spatial dimension, that use a compact stencil (three points in space), see, e.g., [TGB08] for linear and [DFJ03,DFJ04,LK09] for fully nonlinear problems.
More recently, a high-order compact finite difference scheme for (two-dimensional) option pricing models with stochastic volatility has been presented in [DF12a]. This scheme uses a uniform mesh and is fourth order accurate in space and second order accurate in time. Unconditional (von Neumann) stability of the scheme is proved for vanishing correlation. A further study of its stability, indicating unconditional stability also for non-zero correlation, is performed in [DF12b].
In general, the accuracy of a numerical discretisation of (1) for a given number of grid points can be greatly improved by considering a non-uniform mesh. This is particular true for option pricing problems as (1), as typical initial conditions have a discontinuity in their first derivative at S = K, which is the center of the area of interest ('at-the-money').
Our aim in the present paper is to consider extensions of the high-order compact methodology for stochastic volatility models (1) to non-uniform grids. The basic idea of our approach is to introduce a transformation of the partial differential equation from a non-uniform grid to a uniform grid (as, e.g. in [Fou00]). Then, the high-order compact methodology can be applied to this transformed partial differential equation. It turns out, however, that this process is not straightforward as the derivatives of the transformation appear in the truncation error and due to the presence of the cross-derivative terms, one cannot proceed to cancel terms in the truncation error in a similar fashion as in [DF12a] and the derivation of a high-order compact scheme becomes much more involved. Nonetheless, we are able to derive a compact scheme which shows high-order convergence for typical European option pricing problems. Up to the knowledge of the authors, this is the first high-order compact scheme for option pricing in stochastic volatility models on non-uniform grids.
The rest of this paper is organised as follows. In the next section, we transform (1) into a more convenient form. We then derive four new variants of a compact scheme in Section 3. Numerical experiments confirming the high-order convergence for different initial conditions (we consider the case of a European Put option and a European Power Put option) are presented in Section 5. Section 6 concludes.
2 Transformation of the partial differential equation and final condition We focus our attention on the Heston model (1)-(2), although our methodology adapts also to other stochastic volatility models in a natural way (see Remark 2 at the end of Section 3). As usual, we restrict ourselves to the case where the market price of volatility risk λ(S, σ, t) is proportional to σ and choose λ(S, σ, t) = λ 0 σ for some constant λ 0 . This allows to study the problem using the modified parameters which is both convenient and standard practice. For similar reasons, some authors set the market price of volatility risk to zero. The partial differential equation of the Heston-model is then given by imposing an approximative artificial boundary condition at S max . The error caused by approximative boundary conditions imposed on an artificial boundary for a class of Black-Scholes equations has been studied rigorously in [KN00].
The final condition as well as the boundary conditions, which we will discuss separately, depend on the chosen option. In the case of a European Power Put Option we have the final condition For high-order finite difference schemes as proposed in this article, the low regularity of the final condition (6) at the strike S = K may reduce the numerical convergence order in practice. To retain high-order convergence, one can smooth the initial condition carefully (cf. [KTW70]) or shift the numerical grid to avoid the strike falling on a grid point as suggested, for example, in [TR00,DF12a]. In our numerical experiments reported in Section 5 we use the latter approach.
3 Derivation of the high-order compact schemes for the elliptic problem We start by defining a uniform grid in x-and in y-direction, where ∆x = (x max − x min )/N and ∆y = (y max − y min )/M are the step sizes in each direction. With • G we identify the inner points of the grid G. On this grid we denote by U ij the discrete approximation of the continuous solution u in (x i , y j ) ∈ G. Using the standard central difference operator D c x in x-direction and D c y in y-direction, and the standard second-order central difference operator D 2 x in x-direction and D 2 y in y-direction, for k = x, y we have and at the grid points (x i , y j ) for i = 0, . . . , N and j = 0, . . . , M . We call a scheme of high order, if its consistency error is of order O (∆x) 4 for ∆y ∈ O (∆x). If we discretise the higher derivatives u xxxx , u yyyy , u xxxy , u xyyy , u xxx , and u yyy appearing in (9) and (10) with second order accuracy, we obtain a scheme with consistency of order four, since they are all multiplied by factors of order two. If this can be achieved using the compact nine-point computational stencil, the scheme is called high-order compact (HOC).

Auxiliary relations for higher derivatives
We proceed by giving auxiliary relations for the third and fourth order derivatives appearing in (9) and (10). Expressions for the higher derivatives can be obtained by differentiating the partial differential equation (7) in a formal manner without introducing additional error. Differentiating equation (7) with respect to x and then solving for u xxx leads to Using this equation we can calculate a discretisation of A xxx using only points of the nine-point stencil with consistency error of order two using the central difference operators. Differentiating the partial differential equation (7) twice with respect to x and then solving for u xxxx we have The term A xxxx can be discretised at the order two on the compact stencil if equation (11) and the central difference operator are used. Solving equation (12) for u xxxy we obtain In order to find an equation for u yyy we first differentiate the partial differential equation (7) once with respect to y and then solve for u yyy , which leads to The term A yyy can be discretised in a compact manner at the order two using the central difference operators. Differentiating equation (7) twice with respect to y and then solving for u yyyy leads to The term A yyyy can be discretised at the order two on the compact stencil using equation (14) and the central difference operator. Equation (15) is equivalent to Differentiating the partial differential equation (7) once with respect to x and once with respect to y and then solving for u xxxy leads to Using the equations (11) and (14) as well as the central difference operators in x-and y-direction it is possible to discretise A xxxy at the order two on the compact stencil. Solving equation (17) for u xyyy gives Finally, the expression A xyyy can be discretised at the order two on the compact stencil as well.

Derivation of the discrete schemes
In order to derive a discrete scheme we employ equations (9) and (10) in the partial differential equation (7), which gives where is the basis for the derivation of our different discretisation schemes. A 0 is only using the compact stencil.
We have four fourth-order derivatives, namely u xxxx , u yyyy , u xxxy and u xyyy appearing in equation (19), interacting with each other, but only three auxiliary relations to replace these higher derivatives. These relations are given by (12), (15), and (17), which were derived in Section 3.1. This leads to four different versions of the discrete scheme.
In equations (20) to (23) we observe that all these schemes have a formal general consistency error of order two. But on the other hand each version only has one remaining second order term, which is multiplied with either u xxxx , u yyyy , u xxxy , or u xyyy . All other terms are discretised with fourth order accuracy. We call this an essentially high-order compact discretisation. To gauge the overall potential of the four discrete schemes we obtain by neglecting the remaining second-order terms, it is pivotal to understand the behaviour of these terms better. To this end we compute a numerical solution using the (second-order) central difference operator in x-and y-direction directly in equation (7), and obtain by numerical differentiation (approximations of) the higher derivatives u xxxx , u yyyy , u xxxy , and u xyyy appearing in the remaining second order terms.  Figure 1 we see the remainder term for Version 1. This term has by far the highest absolute values. The l 2 -norm of this remainder term is 8.8 × 10 −1 . This indicates that a numerical study of this scheme may not lead to a fourth-order consistency error. On the upper right plot we have the remainder term for Version 2, again without the O((∆x) 2 ) factor. The highest absolute value for this is only about 4 × 10 −3 , so very low when comparing it with the remainder term of Version 1. The l 2 -norm for this plot is 3.1 × 10 −4 , which shows that Version 2 has a significantly higher chance of producing a fourth order consistency error in the numerical study than Version 1. The plot on the lower left side is showing the remainder term of Version 3. This plot has higher values than Version 2, but lower values than Version 1. With a l 2 -norm of 6.6 × 10 −3 it has still a chance to produce a good consistency error. The plot on the lower right shows the remainder term of Version 4. This plot has again very low absolute values which are only up to about 5 × 10 −3 . The l 2 -norm for this remainder term is 3.1 × 10 −4 . This indicates that we have a good chance that Version 4 produces a scheme with fourth-order accuracy.
In the special case that ϕ(x) = x and ∆x = ∆y = h we have (∆y) 2 ≡ (∆x) 2 ϕ 2 x , and all four versions lead to exactly the same HOC scheme, as in this case A xxxy = A xyyy holds. This specific HOC scheme without zoom is discussed in [DF12a].
Remark 2 The derivation of the schemes in this section can be modified to accommodate other stochastic volatility models as, e.g. the GARCH diffusion model (3) or the 3/2-model (4). Using these models the structure of the partial differential equation (1) remains the same, only the coefficients of the derivatives have to be modified accordingly. Similarly, the coefficients of the derivatives in (11)-(18) have to be modified. Substituting these in the modified expression for the truncation error one obtains equivalent approximations as above.
Our conclusion from the results in Figure 1 is that Version 2 and Version 4 seem to be the best choices to obtain small errors. The remainder term for Version 3 still has low values, while Version 1 seems only to be able to produce a second-order scheme. Numerical experiments which we have carried out with all four versions of the scheme indicate that actually Version 3 is leading to the best results in terms of accuracy and stability. Hence, in the remainder of this paper we focus on this particular scheme.
4 High-order compact schemes for the parabolic problem We now consider the parabolic equation (7) with f = −u τ and we denote by U i,j (τ ) the semidiscrete approximation of its solution u(x i , y j , τ ) at time τ .

Semi-discrete schemes
In this section we define the semi-discrete scheme of the form whereK i,j is the coefficient of U i,j (τ ). For the sake of readability we drop the subindex i on the derivatives of ϕ and the subindex j on y, respectively. Analogously we havê as coefficients of ∂ τ U i,j (τ ). With the usage of z ∈ • G we have with n 1 ∈ {i − 1, i, i + 1} and n 2 ∈ {j − 1, j, j + 1}. Thus (24) corresponds to a linear system on • G.

Treatment of the boundary conditions
The first boundary is the boundary x = x min , which corresponds to the boundary at S = 0 of the original problem. For this boundary we have to discount the option price at time T to the appropriate time. Taking into account the transformations τ = T − t and u = e rτ V /K this leads to the Dirichlet boundary condition u(x min , y, τ ) = u(x min , y, 0) for all τ ∈ [0, τ max ] and all y ∈ [y min , y max ].
The next boundary we discuss is the boundary x = x max , which corresponds to the boundary at S = S max of the original problem. For a Power Put with power p we have lim S→∞ V (S, σ, t) = 0, which we approximate at the artificial boundary S max by V S (S max , σ, t) = 0, V SS (S max , σ, t) = 0, V Sσ (S max , σ, t) = 0, V σ (S max , σ, t) = 0 as well as V σσ (S max , σ, t) = 0. Using these approximations in (5) gives V t − rV = 0.
The third boundary to discuss is the boundary y = y min with x / ∈ {x min , x max }, which corresponds to the boundary σ = σ min with S / ∈ {S min , S max }. We will treat this boundary just like the inner of the computational domain, using the equations (25) to (29). This requires the usage of ghost-points U i−1,−1 , U i,−1 and U i+1,−1 when discretising at the points (x i , y 0 ) ∈ G for i = 1, . . . , N − 1. So we need a fourth order accurate expression for the ghost-points U i,−1 for i = 0, . . . , N . We use the following extrapolation formula for i = 0, . . . , N . The same procedure is used for the ghost-points for the matrix M h when using the equations in (30).
The last boundary we discuss is the boundary at boundary y = y max with x / ∈ {x min , x max }, which is corresponding to the boundary σ = σ max with S / ∈ {S min , S max } of the untransformed problem. We treat this boundary similar as the boundary at y min and use equations (25) to (29). The scheme then uses, when discretising at the points (x i , y M ) ∈ G for i = 1, . . . , N − 1, the ghost-points U i−1,M+1 , U i,M+1 and U i+1,M+1 for i = 1, . . . , N − 1. This means that we have to find an expression for the ghost-points U i,M+1 , i = 0, . . . , N . We approximate the values at these ghost-points again using extrapolation, for i = 0, . . . , N . Again, the same procedure is used for the ghost-points for the matrix M h while using the equations in (30).

Time discretization
With the results from the previous sections we obtain a semi-discrete system of the form for each point z of the grid G, which is defined in (8) and ∆x = ∆y = h for some h > 0 is used. The function g(z) has only non-zero values at the boundaries x min and x max .
We use a time grid of the form where the first time steps have step size ∆τ 4 and the following have ∆τ . For these first four time steps, we use the implicit Euler scheme, and obtain with n = 0, 1, 2, 3 for each grid-point z ∈ G. This approach is suggested in [Ran84] when dealing with non-smooth initial conditions. For the following time steps we use a Crank-Nicolson-type time discretisation, leading to with n ≥ 4 on each point z of the grid G. We observe that we have only non-zero values on the compact computational stencil as M x (x) and K x (x) have this property. For the Crank-Nicolson time discretisation this compact scheme has consistency order two in time and four in space for ϕ(x) = x and ρ = 0 or is essentially high-order compact in space otherwise.

Numerical Experiments
In this section we present the results of our numerical experiments for the compact scheme using (25) -(30), whose boundary conditions were derived in Section 4.2. If not stated otherwise, we will use the following default model parameters κ = 1.1, θ = 0.15, v = 0.1, r = ln(1.05), K = 100, T = 0.25.
The initial condition for the European (Power) Put after transformation as in Section 2 is given by where the non-differentiable point of the initial condition is at x K = ϕ −1 (0).

Choice of the zoom function
In our numerical experiments we use the zoom function proposed in [TGB08], with c 1 = asinh(ζŜ min ), c 2 = asinh(ζŜ max ) and ζ > 0. The non-differentiable point of the initial condition hence is at Using the definitions of c 1 and c 2 this can be rearranged tô Hence,Ŝ min can be set by choosing x K in reasonable bounds as well as choosing S max , which giveŝ S max , for a given ζ. The fact that x K can be chosen is very helpful, since if the non-differentiable point is on the grid the numerical convergence order may be reduced to two in practice. Hence, we choose the grid such that the point x K in the middle of two consecutive grid points on the finest grid. This procedure of shifting the grid has been suggested, for example, in [TR00].
In the numerical experiments reported below we choose S min = KeŜ min , S max = 2K, σ min = 0.05, σ max = 0.25. Figure 2 shows the influence of the parameter ζ on the zoom in equation (35), taking into account both transformations,Ŝ = ln (S/K) and x = ϕ −1 (Ŝ). The different values for x K , which depends on ζ, are chosen in such a way that the focus on the values around S = 0 is not too pronounced, compare equation (36). We observe that for smaller values of ζ > 0 there is less zoom. So with ζ → 0 the zoom function is approaching the linear transformation ϕ(x) = (Ŝ max − S min )x +Ŝ min with x ∈ [0, 1]. With a larger value of ζ there is a stronger focus on our area of interest around the exercise price K.
The aim is to find an 'optimal' value for ζ to be used in practical computations. The larger ζ, the smaller the error around K, but on the other hand the error in other parts of the domain increases when having a stronger zoom, because an increasing number of grid points in the area around K automatically results into a decreasing amount of grid points in other areas and vice versa. There has to be a balance between the error in the area around K and the error in other parts of the domain. The overall order of convergence should be looked at to achieve this balance and thus to get a good value for ζ. We expect the numerical convergence order to increase at first with rising ζ and then decrease again after a certain 'optimal' strength of zoom is reached.

Numerical convergence
We now study the numerical errors of the discretisation as h → 0 for fixed parabolic mesh ratio ∆τ /h 2 , using different values for ζ and ρ. We compute an approximation of the solution of the transformed problem, which is given by equation (7), and then transform it back into the original variables. For the relative l 2 -and l ∞ -error plots a reference solution is computed on a fine grid with h ref = 0.003125. For the relative l 2 -errors we use and for the l ∞ -error we use where U ref denotes the reference solution and U is the approximation. We expect the error to behave like O h k for some k. If we plot the logarithm of the error against the logarithm of the number of grid points, the slope of this log-log plot gives the numerical convergence order of the scheme. Due to the initial condition of the transformed problem not being smooth everywhere, we observe that the log-log plots do not always produce a straight line, e.g. for a plain vanilla Put option. For a smooth initial condition the log-log plots of the errors give an almost straight line, e.g. for the Power Put option. The numerical convergence order indicated in the figures below is always computed as the slope of the linear least square fit of the error points. For comparison we additionally plot the results for a standard discretisation (SD), which means that the standard central difference operator is used in (7) as well as In this way all discretisations considered here operate on the same spatial grid and a meaningful comparison can occur. We use ∆τ = 0.4h 2 for all convergence plots, although we note that the dependence of the numerical convergence order on the choice of the parabolic mesh ratio is marginal. This is in line with the results of our numerical stability study reported below in Section 5.3. Looking at the relative l 2 -error we observe that the numerical convergence orders vary from 3.75 to 4.29, which agrees very well with the theoretical order for all zooms. We can also see that the convergence order rises until ζ = 5 and then declines again, so ζ ≈ 5 seems to be the best choice. The lowest relative l 2 -error is always obtained when using ζ = 10. The more useful error in practice is probably the l ∞ -error, as it shows the highest difference between the reference solution and the approximation. When looking at Figure 4 we see that the l ∞ -error and the l 2 -error have a very similar behaviour. The convergence orders vary from 3.00 to 4.20, again having the best order for ζ ≈ 5. When using the finest grid the error for ζ = 5 and ζ = 10 are almost identical, but with rougher grids the error with ζ = 10 is again clearly the lowest. For both error plots we observe that the zoom has its biggest impact when looking at a rough grid, because the error then decreases significantly with an increasing zoom. The HOC discretisations have significantly lower error values and higher convergence orders when comparing them to the standard discretisation. Overall, choosing ζ ≈ 5 for the Heston-Hull-White model (ρ = 0) seems to be the best choice with respect to the convergence order. In Figures 5 and 6 we plot the relative l 2 -and l ∞ -error for a European Put option in the Heston model with ρ = −0.1. This means that the theoretical consistency order is only two, see equation (22). We observe in Figure 5 that the relative l 2 -error varies from 3.40 to 4.14. These values are far above the theoretical consistency order. In fact, using the Version 3 discretisation scheme we obtain a convergence order close to the order using the Heston-Hull-White model. The order of the relative l 2 -error is again rising until ζ = 5 and declining afterwards, but has its lowest values when using ζ = 10. The l ∞ -error in Figure 6 behaves similar to the l ∞ -error in the Heston-Hull-White model. Here the convergence order values vary from 3.00 to 4.09, having its highest value for ζ = 5.
With the finest grid the difference of the error when using ζ = 10 and using ζ = 5 is again very slim. The biggest impact of increasing the zoom in either error plot can be again seen when having a rough grid, because then increasing the zooming leads to significantly lower errors. Similar as in the Heston-Hull-White model the convergence order results are the best when choosing ζ = 5. For both errors we can again see that the essentially high-order compact discretisations have significantly lower error values and higher convergence orders than the standard discretisation.  Figure 7 we can see that the convergence order for the relative l 2 -error varies from 2.92 to 3.84, which is again significantly higher than the theoretical order. The convergence order deteriorates slightly for smaller values of ρ but is still an order better than for the standard discretisation. As expected the best convergence order, which is still very close to four, will be achieved when using ζ = 5. From Figure 8 we find that for the l ∞ -error the convergence order gets lower with lowering the value of ρ. The convergence orders vary from 2.98 to 3.86, where ζ = 5 leads again to the highest value, which is still close to four and thus highly above the theoretical value of the consistency error order. As in the two previous cases the zoom has his highest strengths for the relative l 2 -error as well as for the l ∞ -error when using a very rough grid. For both the relative l 2 -error and the l ∞ -error we can again see that the essentially high-order compact schemes have significantly lower error values and higher convergence orders than the standard discretisation. With the Figures 3 to 8 we recover the numerical observation given in Section 3.2 and can confirm that Version 3 leads to a high-order compact scheme.
For all the discussed European Put options the best results for the convergence order is obtained when using ζ = 5. This value seems to give a good balance between the error around K and the other regions for the zoom. Even though the scheme has a theoretical consistency order equal to four only for the Heston-Hull-White model (ρ = 0), the application showed, that we achieve a numerical convergence order close to four for the Heston model with ρ = 0 as well.
We now consider the case of European Power Put options in the Heston model. The only difference to a plain vanilla European Put is, that the final condition is taken to the power p, see (6), which yields to (34) after transformation. The grid was shifted in a similar manner as above, avoiding x k as a grid point.  In Figures 11 and 12 we can see the convergence orders in the Heston-Hull-White model (ρ = 0) and the Heston model with ρ = −0.4 when p = 3. The differences between the plots are not as big as the theoretical consistency error order may indicate. Even though in the Heston model with ρ = −0.4 the scheme has a theoretical consistency error of order two, it produces a convergence order from 3.50 to 3.69 depending on the zoom strength ζ, whereas the orders in the Heston-Hull-White model with ρ = 0, where we have a theoretical consistency order of four, vary from 4.04 to 4.10. In both situations the standard discretisation is outperformed on behalf of convergence order and error values.

Numerical stability study
In the particular case of a uniform grid, i.e. ϕ(x) = x, the scheme developed here reduces to the high-order compact scheme presented in [DF12a], where unconditional (von Neumann) stability is proved for ρ = 0. An additional stability analysis performed in [DF12b] suggests that the scheme is also unconditionally stable for general choice of parameters. For the present scheme on a non-uniform grid, a similar von Neumann analysis, analytical or numerical, appears to be out of reach as the expression for the amplification factor is formidable and consists of high-order polynomials in a two-digit number of variables. To validate the stability of the scheme for general parameters, we therefore perform additional numerical stability tests. We remark that in our numerical experiments we observe a stable behaviour throughout. We compute numerical solutions for varying values of the parabolic mesh ratio c = ∆τ /h 2 and the mesh width h. Plotting the associated relative l 2 -norm errors in the plane should allow us to detect stability restrictions depending on c or oscillations that occur for high cell Reynolds number (large h). This approach for a numerical stability study was also used in [DF12a,DFJ03].
We show results for the European Put option in the Heston Model only, since the Power Puts only differ in the initial conditions and give similar results. For our stability plots we use c = k/10 with k = 1, . . . , 10, and a descending sequence of grid points in x-direction, starting with six grid points (since x ∈ [0, 1] it follows h ≤ 0.2), and doubling the number of points (halving h) in each step. The zoom parameter ζ = 5 is used.

Conclusion
We have presented new high-order compact finite difference schemes for option pricing under stochastic volatility on non-uniform grids. The resulting schemes are fourth-order accurate in space and second-order accurate in time for vanishing correlation. In our numerical convergence study we obtain high-order numerical convergence also for non-zero correlation and non-smooth payoffs which are typical in option pricing. In all numerical experiments a comparative standard second-order discretisation is significantly outperformed. We have conducted a numerical stability study which seems to indicate unconditional stability of the scheme. In our numerical experiments we observe a stable behaviour for all choices of parameters.
It would be interesting to consider extensions of this scheme to the American option pricing problem, where early exercise of the option is possible. In this case, one has to solve a free boundary problem. It can be written as a linear complementarity problem which could be discretised using the schemes given here. To retain the high-order convergence one would need to combine the highorder discretisation with a high-order resolution of the free boundary. This extension is beyond the scope of the present paper, and we leave it for future research.

A Coefficients for Version 2 and Version 4
In this section we give the coefficients of the semi-discrete schemes for Version 2 and Version 4. We do not include the coefficients for Version 1 as this version always resulted into a second-order numerical convergence error in the numerical study.

A.1 Coefficients for Version 2
When discretising equation (21) with the central difference operator in x-and in y-direction, we get the following coefficients for the Version 2 schemê