All articles published by MDPI are made immediately available worldwide under an open access license. No special

permission is required to reuse all or part of the article published by MDPI, including figures and tables. For

articles published under an open access Creative Common CC BY license, any part of the article may be reused without

permission provided that the original article is clearly cited. For more information, please refer to

https://www.mdpi.com/openaccess.

Feature papers represent the most advanced research with significant potential for high impact in the field. A Feature

Paper should be a substantial original Article that involves several techniques or approaches, provides an outlook for

future research directions and describes possible research applications.

Feature papers are submitted upon individual invitation or recommendation by the scientific editors and must receive

positive feedback from the reviewers.

Editor’s Choice articles are based on recommendations by the scientific editors of MDPI journals from around the world.

Editors select a small number of articles recently published in the journal that they believe will be particularly

interesting to readers, or important in the respective research area. The aim is to provide a snapshot of some of the

most exciting work published in the various research areas of the journal.

Inspired by the article Weak Convergence Rate of a Time-Discrete Scheme for the Heston Stochastic Volatility Model, Chao Zheng, SIAM Journal on Numerical Analysis 2017, 55:3, 1243–1263, we studied the weak error of discretization schemes for the Heston model, which are based on exact simulation of the underlying volatility process. Both for an Euler- and a trapezoidal-type scheme for the log-asset price, we established weak order one for smooth payoffs without any assumptions on the Feller index of the volatility process. In our analysis, we also observed the usual trade off between the smoothness assumption on the payoff and the restriction on the Feller index. Moreover, we provided error expansions, which could be used to construct second order schemes via extrapolation. In this paper, we illustrate our theoretical findings by several numerical examples.

The Heston Model Heston (1993) is a widely used stochastic volatility model to price financial options. It consists of two stochastic differential equations (SDEs) for an asset price process S and its volatility V:

with , , , and independent Brownian motions , , which are defined on a filtered probability space , where the filtration satisfies the usual conditions. It is a simple and popular extension of the Black–Scholes model where the volatility of the asset was assumed to be constant. As a consequence, the Heston Model takes the asymmetry and excess kurtosis of financial asset returns into account which are typically observed in real market data. The volatility is given by the so-called Cox–Ingersoll–Ross process (CIR). Its Feller index will be an important parameter for our results. Throughout this article, the initial values , are assumed to be deterministic.

To price options with maturity at time T, one is interested in the value of

where is the payoff function. Closed formulae for are rarely known and often Monte Carlo methods are applied, for which in turn the simulation of is required. Usually, the log-Heston model instead of the Heston model is considered in numerical practice. This yields the SDE

and the exponential is then incorporated in the payoff, i.e., g is replaced by with .

A second order discretization scheme for the log-Heston model has been introduced in Andersen (2008) and analyzed in Zheng (2017). The so-called Broadie-Kaya trick and a removal of the drift, detailed in Section 3.1, reduce the simulation of the log-Heston model to the joint simulation of

Moreover, since the transition density of the CIR process follows a non-central chi-square distribution, it can be simulated exactly. Trapezoidal discretizations of the first component lead to the trapezoidal scheme

where , and . This discretization avoids in particular the cumbersome exact simulation of the integrated volatility. Zheng (2017) establishes weak order two for polynomial test functions by transferring the error analysis to that of a trapezoidal rule for multidimensional deterministic integrals. Our original intention was to extend this result to a larger class of test functions f by using the Kolmogorov PDE approach. However, the required Itō-Taylor expansions turned out to be not feasible. So, instead, we analyzed the following two semi-exact discretization schemes: the Euler-type scheme

and the semi-trapezoidal scheme

In both schemes, the CIR process is simulated exactly. In our opinion, the analysis of these schemes gives valuable insights in the weak error analysis of discretization schemes for the log-Heston model and is also a good starting point for the analysis of full Euler-type discretization schemes.

Our error analysis relies on two regularity results for the Heston PDE (Briani et al. (2018); Feehan and Pop (2013)), the Kolmogorov PDE approach for the weak error analysis from Talay and Tubaro (1990), and Malliavin calculus. We also observe the usual trade off between the smoothness assumption on the payoff and the restriction on the Feller index. For payoffs of lower smoothness, a restriction on the Feller index is required, which arises from the use of Malliavin calculus tools.

In the following, we use the notation

for the maximal step size and the usual notations for the spaces of differentiable functions. In particular, the subscript c denotes compact support and denotes polynomial growth. In addition, see Section 3.1. The results of Feehan and Pop (2013) require compact support of the test functions f, while the results of Briani et al. (2018) allow polynomial growth but require higher smoothness for f.

Theorem1.

Let . (i) If and , then both schemes satisfy

(ii) If , then both schemes satisfy

Assuming more smoothness of f, we obtain more detailed results:

Theorem2.

Suppose that . (i) Then, the Euler scheme (4) satisfies

where

and

for.

In particular, for an equidistant discretization with, , we have

Here, u denotes the solution of the associated Kolmogorov PDE; see Equation (7).

In particular, for an equidistant discretization , , it holds

Here, u denotes again the solution of the associated Kolmogorov PDE; see Equation (7).

Thus, the semi-trapezoidal rule eliminates the first two terms of the error expansion of the Euler scheme.

Remarks

Remark1.

We expect that the error expansions for an equidistant discretization for both schemes satisfy

However, to establish this, we would require error estimates for functionals of the type with , which are uniform in λ. (Compare, e.g., Proposition 2 in Talay and Tubaro (1990).) This, in turn, would require uniform regularity estimates for the Heston PDE, which are not available at the moment.

Remark2.

Property (6) allows to construct a second order scheme via extrapolation: If (6) holds, then

satisfies

where uses the stepsize and the stepsize .

Remark3.

We require smoothness assumptions for f that are not met by the payoffs in practice, which are at most Lipschitz continuous or even discontinuous. However, this is a typical problem for weak approximation of SDEs as the Heston SDE, which do not satisfy the so-called standard assumptions on the coefficients. In Bally and Talay (1996), only bounded and measurable test functions f are treated assuming uniform hypoellipticity of the coefficients of the SDE. However, the Heston model does not satisfy this property. An adaptation of the strategy of Bally and Talay (1996) to the Heston model yields strong assumption on the Feller index (see Altmayer (2015)), which we want to avoid here.

Remark4.

Schemes built on the Broadie-Kaya trick, i.e., Equation (3), have a different structure than schemes which arise by a direct discretization of the log-Heston model as, e.g., the schemes studied in Altmayer and Neuenkirch (2017); Lord et al. (2009). For example, the so-called absorbed Euler discretization reads as

Here, the volatility is discretized by an Euler scheme, a fix for retaining the positivity is introduced by using the positive part, and the equation for the log-Heston price is discretized instead of the one for .

Remark5.

The Broadie-Kaya trick is a particular case of a more general transformation procedure, which has been introduced in Cui et al. (2018) for a general class of stochastic volatility models. In addition, in Cui et al. (2018), the weak convergence of a Markov chain approximation for these equations is established, which had been introduced in Cui et al. (2020). Markov chain approximations have been also studied in Briani et al. (2018) for the Heston and Bates model and are an alternative to a classical discretization of stochastic differential equations. In particular, for pricing American options, they can be beneficial.

2. Numerical Results

In this section, we will test numerically whether the convergence rates for the Euler Scheme (4) and the Semi-Trapezoidal Scheme (5) are attained even under milder assumptions than those from Theorems 1 and 2. We use the following model parameters:

Model 1: ;

Model 2: ;

Model 3: .

The Feller index is in Model 1, in Model 2, and in Model 3. For each model, we use the following payoff functions:

1.

European Call: ;

2.

European Put: ;

3.

Indicator: .

Note that none of these payoffs satisfies the assumptions of our Theorems. Thus, the presented numerical experiments explore whether the Theorems are valid under milder assumptions. In order to measure the weak error rate, we simulated independent copies , , of to estimate

by

for each combination of model parameters, functional and number of steps where . The number of Monte Carlo samples is chosen in such a way that the Monte Carlo error is sufficiently small enough, i.e., does not dominate the theoretically expected convergence rates. The Monte Carlo mean of these samples was then compared to a reference solution , i.e.,

and the error is plotted in Figure 1, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15, Figure 16, Figure 17 and Figure 18. We then measured the rate of convergence, i.e., the decay rate of , by the slope of a least-squares fit in logarithmic coordinates. The reference solutions can be computed with sufficiently high accuracy from semi-explicit formulae via Fourier methods. In particular, the put price can be calculated from the call-price formula given in Heston (1993) via the put-call-parity. The price of the digital option can be computed from the probability given in Heston (1993); it equals . Additionally to the Euler and Semi-Trapezoidal scheme, we simulated the Trapezoidal scheme as in Zheng (2017) and the two extrapolation schemes from Remark 2. Moreover, to present a broader picture we estimated the weak error order of two Euler-type discretizations of the full Heston Model, the Full Truncation Euler (FTE) as in Lord et al. (2009), and the Symmetrized Euler as in Bossy and Diop (2015). To clarify things, we show two plots for each combination of model parameters and functional: one with the suspected order one schemes (Euler, Semi-Trapezoidal, FTE, and Symmetrized Euler) and one with the suspected order two schemes (Trapezoidal, Extrapolated Euler, and Extrapolated Semi-Trapezoidal).

All “Order 1” schemes seem to have a very regular convergence behavior except for the Semi-Trapezoidal scheme for the Indicator, which could be explained by the low absolute error. Especially for the Call and the Indicator, both schemes from Theorem 1 seem to have very high weak convergence rates. Because of the Feller index of 0.63 in this model, this indicates that the assertion of Theorems 1 and 2 could hold under weaker assumptions. The extremely low estimated convergence rate for the Semi-Trapezoidal scheme in combination with the Put could be due to the low error. The estimated weak error order of the FTE scheme is noticeably higher than 1, whereas the Symmetrized Euler has low convergence rates. The convergence behavior of the “Order 2” schemes is a bit less regular. The Extrapolated Euler scheme seems to converge with order 2 for all payoff functions, whereas the Extrapolated Semi-Trapezoidal scheme seem to have only order 1 for the Indicator. But, again, we notice that the error for just 2 discretization steps already starts at around 2−10, which is extremely low.

2.2. Model 2

Here, we have an even lower Feller index of . We can see that the estimated convergence rates for all “Order 1” schemes are lower than before, see Table 2. However, the Semi-Trapezoidal scheme and the FTE scheme seem to converge with order 1. The convergence behavior is still quite regular as we can see in Figure 7, Figure 9, and Figure 11. In absolute terms, the errors of the schemes from Theorem 1 are the lowest, especially for Put and Indicator. Looking at the “Order 2” schemes, the Trapezoidal discretization still shows an estimated weak convergence rate of around 2, whereas the two extrapolation schemes show a weaker performance. But, especially for the Indicator, all three schemes seem to have a very low error and a quite regular convergence behavior.

2.3. Model 3

Here, we have the highest Feller Index with . It is, therefore, a bit surprising that the Euler scheme seems to have a convergence rate of less than 1 in this case. In general, the errors for the “Order 1” schemes show a more irregular behavior, as can be seen from Figure 13, Figure 15, and Figure 17. The Semi-Trapezoidal and the FTE scheme work especially well in this scenario as we can see in Table 3. This is also the only case where the Symmetrized Euler shows an estimated convergence order of around 1. The extrapolation definitely improves the convergence rate of the Euler scheme with order 2 for the Indicator, but this is not the case for the Semi-Trapezoidal scheme.

2.4. Computational Times

The computational times show the expected behavior, i.e., the simulation times for the semi-exact schemes increase as the Feller index decreases. See Table 4 and Table 5. This is a well known feature of the MATLAB-generator ncx2rnd for the non-central chi-square distribution, which we used. (All simulations were carried out in MATLAB.)

2.5. Conclusions

Except for the Euler scheme for the Call in Model 3, the simulation studies support the conjecture that the convergence rates of Theorems 1 and 2 hold under weaker assumptions. For the mentioned behavior of the Euler scheme, we do not have an explanation, except the possibly pre-asymptotic step sizes. For the extrapolated schemes, which might have order two, the situation is less clear. Since the behavior of the trapezoidal scheme is regular, a too large Monte Carlo error seems an unlikely explanation. Explanations could be again the pre-asymptotic step sizes or, in fact, the non-smoothness of the considered payoffs.

3. Auxiliary Results

In this section, we will collect and establish, respectively, several auxiliary results for the weak error analysis.

3.1. Kolmogorov PDE

Recall that the stochastic integral equations for the log-Heston model for read as

Now, we apply the so-called Broadie-Kaya trick from Broadie and Kaya (2006). We can rearrange the second equation:

Then, we plug this equation into the first one:

Without loss of generality, we can neglect the non-integral part in , since we have

with given below. To get the Kolmogorov backward PDE, we look at the following integral equations:

We set

and obtain for bounded and continuous the Kolmogorov backward PDE by an application of the Feynman-Kac Theorem (see, e.g., Theorem 5.7.6 in Karatzas and Shreve (1991)):

In our error analysis, we will follow the now classical approach of Talay and Tubaro (1990), which exploits the regularity of the Kolmogorov backward PDE. For the latter we will rely on the works of Feehan and Pop (2013) and Briani et al. (2018). To state these regularity results, we will need the following notation:

For a multi-index , we define and for , we define . Moreover, we denote by the standard Euclidean norm in . Let be a domain and . The set is the set of all real-valued functions on which are q-times continuously differentiable. For , we denote by the set of all functions from in which partial derivatives of order q are Hölder-continuous of order , and is the set of all functions from , who have compact support. Moreover, is the set of functions such that there exist for which

Finally, we denote by the set of functions such that there exist for which

The work of Feehan and Pop deals with general degenerated parabolic equations and establishes a-priori regularity estimates for them. In the context of Equation (7), the main result of Feehan and Pop (2013), i.e., Theorem 1.1, reads as follows:

Theorem3.

Let and . Then, there exists a constant , depending only on and σ such that the solution u of PDE (7) satisfies

So, under the above assumptions on f, the solution u and the first order derivatives are bounded. Moreover, the second order derivatives are also bounded, if they are damped by v for .

Assuming more smoothness on f, we can achieve more regularity for u using the above result, at least for the partial derivatives with respect to x. Set

This is well defined: by continuity and boundedness of and dominated convergence we have

with

An analogous calculation for shows that . Thus, is also bounded, if . Moreover, fulfills the Kolmogorov backward PDE

while fulfills the same PDE with terminal condition

Applying Theorem 3 now to and , we obtain the following additional bounds (case (ii)) for the derivatives of u:

Corollary1.

(i) Let and . Then, there exists a constant , depending only on and σ such that the solution u of PDE (7) satisfies

(ii) Let and . Then, there exists a constant , depending only on and σ such that the solution u of PDE (7) satisfies

The recent work of Briani et al. is a specialized approach for the log-Bates model, of which the log-Heston model is a particular case. In our setting, they obtain in Proposition 5.3 and Remark 5.4 of Briani et al. (2018) the following:

Theorem4.

Let , and suppose that . Then, the solution u of PDE (7) satisfies .

In contrast to the results of Feehan and Pop, the result of Briani et al. requires more smoothness of f but allows polynomial growth instead of compact support.

3.2. Properties of the CIR Process

We recall here the following estimates for the CIR process, which are well known or can be found in Hurd and Kuznetsov (2008).

Lemma1.

(1) We have

for all and

(2) For all , there exist constants , depending only on and , such that

We will need the following bound on the growth of the -norm of a specific stochastic integral of the CIR process:

Lemma2.

For all , it holds that

Proof.

With the Burkholder-Davis-Gundy inequality and the Hölder inequality, we have

for all . The assertion now follows from Lemma 1 (1). □

3.3. Malliavin Calculus

When working with low smoothness assumptions on f, we will use a Malliavin integration by parts procedure to establish weak convergence order one. As in Altmayer and Neuenkirch (2017), this paragraph gives a short introduction into Malliavin calculus; for more details, we refer to Nualart (1995).

Malliavin calculus adds a derivative operator to stochastic analysis. Basically, if Y is a random variable and a two-dimensional Brownian motion, then the Malliavin derivative measures the dependence of Y on . The Malliavin derivative is defined by a standard extension procedure: Let be the set of smooth random variables of the form

with bounded with bounded derivatives, , , and the stochastic integrals

The derivative operator D of such a smooth random variable is defined as

This operator is closable from into with and the Sobolev space denotes the closure of with respect to the norm

In particular, if denotes the first component of the Malliavin derivative, i.e., the derivative with respect to W, we have

and vice versa for the derivative with respect to B, i.e.,

This, in particular, implies that, if is independent of B, then .

For the CIR process, we will, therefore, have that for all .

The derivative operator follows rules similar to ordinary calculus.

Proposition1.

Let be a random variable with components in . If

(i)

is in ,

(ii)

,

(iii)

,

then the chain rule holds: and

For example, for a random variable and with bounded derivative, the chain rule reads as

Another simple example for the application of this chain rule is

The divergence operator is the adjoint of the derivative operator. If a random variable belongs to , the domain of the divergence operator, then is defined by the duality—also called integration by parts—relationship

If u is adapted to the canonical filtration generated by and satisfies , then and coincides with the Itō integral . For the Malliavin regularity of the CIR process, the following is well known. See, e.g., Proposition 4.5 and Theorem 4.6 in Altmayer (2015) or Proposition 4.1 in Alos and Ewald (2008).

under the assumption with differentiable and bounded, see Proposition 4.1 in Altmayer and Neuenkirch (2015). Indeed, using and , and the chain rule, i.e.,

we have

where the first equality is due to the integration by parts formula.

In Lemmas 5 and 9, we will establish discrete counterparts for this integration by parts result, i.e., on the level of the approximation schemes. In this context, we will also need the Malliavin differentiability of . Since

we obtain

by exchanging the Riemann integral and the Malliavin derivative (via a standard approximation argument for the Riemann integral, Lemma 3 and Lemma 1.2.3 in Nualart (1995)) and the independence of and B. Thus, we can conclude that

3.4. Properties of the Euler Discretization

Recall that the Euler discretization of the price process is given by

with . We extend this discretization in every interval as the following Itō process:

Here, we have set , and .

We have the following result on the Malliavin regularity of the Euler discretization:

Lemma4.

Let and . Then, , and we have

Proof.

We have

and

Following the steps of the proof of Lemma 3.5 from Altmayer and Neuenkirch Altmayer and Neuenkirch (2017), we then have exploiting that and under the assumption . The chain rule from Proposition 1 yields

and

□

Note that we write, in the following, instead of to unify the notation. With the above result, we can express without the second order derivative of u, which will be needed later on.

Lemma5.

Let . Under the assumptions of Theorem 3 and , we have

Proof.

To avoid stronger restrictions on the Feller index we will use a localization procedure. So, for , let be a function such that

1.

is continuously differentiable with bounded derivative,

2.

on ,

3.

on ,

4.

on .

Since and B are independent, the chain rule from Proposition 1 implies

with . Recall the integration by parts formula from Equation (8), i.e.,

where we now choose

Before we can apply the integration by parts rule, we need to check whether

for . We deduced these terms by using again the chain rule for . Note that the properties of the localizing function and Theorem 3 imply that

are all uniformly bounded in . So, Equation (10) holds, then, due to Lemma 1, Lemma 3, Equation (9), and Lemma 4.

Since is also well-defined by Lemma 1 due to , we obtain now

Due to Corollary 1 (i), not only but also is bounded. Since almost surely for and for all , the assertion follows now by dominated convergence using the Itô-isometry and again Lemma 1. □

We also need the following -convergence result:

Lemma6.

Let . There exists a constant , depending only on and , such that

Proof.

We have

Assume without loss of generality that . Jensen’s inequality and the Burkholder-Davis-Gundy inequality now imply that there exists a constant , depending only on p, T, the parameters of the CIR process, and , such that

Since for , the assertion follows from Lemma 1. □

Straightforward calculations also yield the following -smoothness result for the Euler-type scheme:

Lemma7.

Let . There exists a constant , depending only on , and , such that

for all .

3.5. Properties of the Semi-Trapezoidal Rule

Recall that our semi-trapezoidal rule reads as

Again, we write the scheme as a time-continuous process:

Expanding the last term with Itō’s lemma, we obtain

with

Here, we have set again , , , and we also write again , instead of , to unify the notation.

We need the following result on the Malliavin regularity of the semi-trapezoidal scheme:

Lemma8.

Let and . Then, we have and

Proof.

We already know that and . We can write as

with

Following the steps of the proof of Lemma 3.5 from Altmayer and Neuenkirch (2017), we then also have . The chain rule from Proposition 1 yields

and

□

Note that the partial Malliavin derivative with respect to B for the Euler and the semi-trapezoidal scheme coincide. So, by analogous calculations as for the Euler scheme, we obtain the following integration by parts result:

Lemma9.

Let . Under the assumptions of Theorem 3 and , we have

By similar calculations as for the Euler scheme, we also have:

Lemma10.

Let . There exists a constant , depending only on , and , such that

Lemma11.

Let . There exists a constant , depending only on , and , such that

for all .

4. Proof of Theorem 1

We address both schemes and the different assumptions in separate subsections. Constants, which are, in particular, independent of the maximal stepsize

and depend only , and , will be denoted by c, regardless of their value.

4.1. The Euler Scheme: Expanding the Error

Since and , the weak error is a telescoping sum of local errors:

With the Itō formula and the Kolmogorov backward PDE evaluated at , we obtain

Since

we have with

By Theorem 3 and Corollary 1, we have that and are bounded. So, Lemma 1 implies that

and

Moreover, with the law of total expectation and the adaptedness of and , we have

due to the martingale property of the Itō integral. Therefore, we can write

and obtain

In the same way, we have

Summarizing this preliminary part, we have obtained

where

4.2. The Euler Scheme: Case (i)

So, it remains to analyze and under the regularity of Theorem 3 (i). We start with . The mean value theorem and

and give

where we used

With the Minkowski inequality and Lemma 1, it holds that

for , where c is in particular independent of . Finally, with the Minkowski inequality, the Burkholder-Davis-Gundy inequality, Lemma 1, and Lemma 7, we obtain for all that

where c is in particular independent of . The Hölder inequality then gives

For , we will use the integration by parts rule to get rid of the second order derivative. Otherwise, direct estimation would only lead to weak order . First, recall that, by Lemma 5, we have

Moreover, note that we also have

and recall that

Thus, we can write

with .

Before we analyze this expression further, it remains to show (14). Using the law of total expectation, the adaptedness of , and of the Itō integrals, we have

Since

due to the properties of the Itō integral, Equation (14) follows.

Using the mean value theorem in Equation (15), we obtain

Therefore,

with

By Lemmas 1 and 2, it holds that

for and . So, the Hölder inequality leads to

with and c in particular independent of .

With the Cauchy-Schwarz, Burkholder-Davis-Gundy, and Minkowski inequalities for , it follows that

With the Hölder inequality, we now have

Therefore,

and, since is Riemann-integrable, we obtain

which concludes the proof of this part.

4.3. The Euler Scheme: Case (ii)

Starting from Equation (11) and using now the bounds of Corollary 1 for , , , and , the assertion follows from a direct application of the mean value theorem to (12) and (13), together with the Lemmata 1 and 7.

4.4. Semi-Trapezoidal Rule: Expanding the Error

We look again at the telescoping sum of local errors

Recall that

with

With the Itō formula and the Kolmogorov backward PDE evaluated at , we have

Using again

we obtain

with

Since by Theorem 3, we have

using Lemma 1. Moreover, since and are bounded by Theorem 3 and Corollary 1 (i), we obtain similar to the calculations for the Euler scheme that

and

with

4.5. Semi-Trapezoidal Rule: Case (i)

Since Lemma 8 gives and

we can proceed here in the same way as for the Euler scheme by using the Lemmata 9 and 11.

4.6. Semi-Trapezoidal Rule: Case (ii)

Starting from (16), the assertion of this case follows from a direct application of the mean value theorem to (17) and (18) using the regularity results from Corollary 1, together with the Lemmata 1 and 11.

5. Proof of Theorem 2

Now, we derive the error expansion under the regularity of Theorem 4 with , i.e., we have .

5.1. Euler Scheme: Preliminaries

By the Lemmata 1, 4, and 7, we have that

and

for all . Using the Burkholder-Davis-Gundy, Hölder, and Minkowski inequalities, we also have

and

for all . We will use this in the following at several places without explicitly mentioning it.

If higher derivatives of u are available, then we can analyze

and

via another application of Itō’s lemma. So, let be a -function that fulfills the backward PDE (7). In particular, the partial derivatives of u up to order two are such functions. Itō’s formula and the Kolmogorov backward PDE (7) now give

If and have polynomial growth, then an application of the Itō isometry and the martingale property of the Itō integral yield

where

If and have polynomial growth, then an application of Hölder’s inequality and the Itō isometry yield

and so it follows

Since we have

again, by the properties of the Itō integral, we finally obtain by Hölder’s inequality and the Burkholder-Davis-Gundy inequality that

Thus, we can conclude that

for and , if the derivatives up to order four of u have polynomial growth.

5.2. Euler Scheme: Conclusion

Setting now in Equation (21), we have from (19) that

Replacing now the function k by in Equation (21), we arrive from (20) at

Summarizing, we have shown that

and

where

An application of the mean value theorem, the polynomial growth of the derivatives of u, the Minkowski inequality, the Hölder inequality, and the Lemmata 1, 6 yields for that

Note here that implies that and are well-defined, have polynomial growth, and are continuous.

Thus, for an equidistant discretization , , we have

Since

for , this concludes the proof of Theorem 2 (i).

5.3. Semi-Trapezoidal Scheme: Preliminaries

By the Lemmata 1, 8, and 11, we have that

and

for all . Using the Burkholder-Davis-Gundy, Hölder, and Minkowski inequalities, we also have

and

for all .

We will use this in the following at several places without explicitly mentioning it.

We will now take also a closer look at the error of the semi-trapezoidal discretization for . Recall that

We can again use the Itō formula and the Kolmogorov backward PDE (7) evaluated at and obtain for a -function k, which fulfills the PDE (7), that

with

Analogous calculations as for the Euler scheme yield that

for and under the assumption .

5.4. Semi-Trapezoidal Rule: Calculations for , , and

Rewriting the terms of using (23) for the last term gives

Applying, again, (23) with s instead of as the lower bound of the integral to the second summand of the first term, using the polynomial growth of the derivatives of u and Hölder’s inequality, we also have

and so

Adding yields

However, Itō’s formula gives for sufficiently smooth that

with

where

and

Since , we can apply this to and taking expectations gives then

So, we end up with

Looking at , the last term is already of third order:

An application of the mean value theorem, the polynomial growth of the derivatives of u, the Minkowski inequality, the Hölder inequality, and the Lemmata 1, 10 yields for that

In particular, for an equidistant discretization , , we have

and the convergence

for concludes the proof of Theorem 2 (ii).

Author Contributions

All authors have contributed substantially to this manuscript. All authors have read and agreed to the published version of the manuscript.

Funding

DFG Research Training Group 1953 “Statistical Modeling of Complex Systems”.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data sharing not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

Alos, Elisa, and Christian-Oliver Ewald. 2008. Malliavin differentiability of the Heston volatility and applications to option pricing. Advances in Applied Probability 40: 144–62. [Google Scholar] [CrossRef]

Altmayer, Martin. 2015. Quadrature of Discontinuous SDE Functionals Using Malliavin Integration by Parts. Ph.D. dissertation, University of Mannheim, Germany. [Google Scholar]

Altmayer, Martin, and Andreas Neuenkirch. 2015. Multilevel Monte Carlo quadrature of discontinuous payoffs in the generalized Heston model using Malliavin integration by parts. SIAM Journal on Financial Mathematics 6: 22–52. [Google Scholar]

Altmayer, Martin, and Andreas Neuenkirch. 2017. Discretising the Heston model: An analysis of the weak convergence rate. IMA Journal of Numerical Analysis 37: 1930–60. [Google Scholar] [CrossRef]

Andersen, Leif. 2008. Simple and efficient simulation of the Heston stochastic volatility model. Journal of Computational Finance 11: 29–50. [Google Scholar]

Bally, Vlad, and Denis Talay. 1996. The law of the Euler scheme for stochastic differential equations. I: Convergence rate of the distribution function. Probability Theory and Related Fields 104: 43–60. [Google Scholar] [CrossRef]

Bossy, Mireille, and Awa Diop. 2015. Weak convergence analysis of the symmetrized Euler scheme for one dimensional SDEs with diffusion coefficient |x|a, a ∈ [1/2,1). arXiv, arXiv:1508.04573. [Google Scholar]

Briani, Maya, Lucia Caramellino, and Giulia Terenzi. 2018. Convergence rate of Markov chains and hybrid numerical schemes to jump-diffusions with application to the Bates model. arXiv, arXiv:1809.10545. [Google Scholar]

Broadie, Mark, and Özgür Kaya. 2006. Exact simulation of stochastic volatility and other affine jump diffusion processes. Operations Research 54: 217–31. [Google Scholar] [CrossRef]

Coskun, Sema, and Ralf Korn. 2018. Pricing barrier options in the Heston model using the Heath–Platen estimator. Monte Carlo Methods and Applications 24: 29–41. [Google Scholar]

Cui, Zhenyu, Justin Kirkby, and Duy Nguyen. 2018. A general valuation framework for SABR and stochastic local volatility models. SIAM Journal on Financial Mathematics 9: 520–63. [Google Scholar]

Cui, Zhenyu, Justin Kirkby, and Duy Nguyen. 2020. Efficient simulation of generalized SABR and stochastic local volatility models based on markov chain approximations. European Journal of Operational Research. [Google Scholar] [CrossRef]

Feehan, Paul, and Camelia Pop. 2013. A Schauder approach to degenerate-parabolic partial differential equations with unbounded coefficients. Journal of Differential Equations 254: 4401–45. [Google Scholar] [CrossRef]

Glasserman, Paul, and Kyoung-Kuk Kim. 2011. Gamma expansion of the Heston stochastic volatility model. Finance and Stochastics 15: 267–96. [Google Scholar] [CrossRef][Green Version]

Heston, Steven. 1993. A closed-form solution for options with stochastic volatility with applications to bond and currency options. The Review of Financial Studies 6: 327–43. [Google Scholar] [CrossRef][Green Version]

Hurd, Thomas, and Alexey Kuznetsov. 2008. Explicit formulas for Laplace transforms of stochastic integrals. Markov Processes and Related Fields 14: 277–90. [Google Scholar]

Kahl, Christian, and Peter Jäckel. 2006. Fast strong approximation Monte-Carlo schemes for stochastic volatility models. Quantitative Finance 6: 513–36. [Google Scholar]

Karatzas, Ioannis, and Steven Shreve. 1991. Brownian Motion and Stochastic Calculus, 2nd ed. New York: Springer. [Google Scholar]

Malham, Simon, and Anke Wiese. 2013. Chi-square simulation of the CIR process and the Heston model. International Journal of Theoretical and Applied Finance 16: 1350014. [Google Scholar]

Nualart, David. 1995. The Malliavin Calculus and Related Topics. New York: Springer. [Google Scholar]

Lord, Roger, Remmert Koekkoek, and Dick van Dijk. 2009. A comparison of biased simulation schemes for stochastic volatility models. Quantitative Finance 10: 177–94. [Google Scholar]

Smith, Robert. 2007. An almost exact simulation method for the Heston model. Journal of Computational Finance 11: 115–25. [Google Scholar]

Talay, Denis, and Luciano Tubaro. 1990. Expansion of the global error for numerical schemes solving stochastic differential equations. Stochastic Analysis and Applications 8: 483–509. [Google Scholar]

Zheng, Chao. 2017. Weak convergence rate of a time-discrete scheme for the Heston stochastic volatility. SIAM Journal on Numerical Analysis 55: 1243–63. [Google Scholar] [CrossRef]

Figure 1.

Call Model 1.

Figure 1.

Call Model 1.

Figure 2.

Call Model 1.

Figure 2.

Call Model 1.

Figure 3.

Put Model 1.

Figure 3.

Put Model 1.

Figure 4.

Put Model 1.

Figure 4.

Put Model 1.

Figure 5.

Indicator Model 1.

Figure 5.

Indicator Model 1.

Figure 6.

Indicator Model 1.

Figure 6.

Indicator Model 1.

Figure 7.

Call Model 2.

Figure 7.

Call Model 2.

Figure 8.

Call Model 2.

Figure 8.

Call Model 2.

Figure 9.

Put Model 2.

Figure 9.

Put Model 2.

Figure 10.

Put Model 2.

Figure 10.

Put Model 2.

Figure 11.

Indicator Model 2.

Figure 11.

Indicator Model 2.

Figure 12.

Indicator Model 2.

Figure 12.

Indicator Model 2.

Figure 13.

Call Model 3.

Figure 13.

Call Model 3.

Figure 14.

Call Model 3.

Figure 14.

Call Model 3.

Figure 15.

Put Model 3.

Figure 15.

Put Model 3.

Figure 16.

Put Model 3.

Figure 16.

Put Model 3.

Figure 17.

Indicator Model 3.

Figure 17.

Indicator Model 3.

Figure 18.

Indicator Model 3.

Figure 18.

Indicator Model 3.

Table 1.

Measured convergence rates Model 1.

Table 1.

Measured convergence rates Model 1.

Method

Call

Put

Indicator

Euler

1.5252

0.9492

1.1870

Semi-Trapezoidal

2.0174

0.2857

1.8343

FTE

1.5205

1.5205

1.2847

Symmetrized Euler

0.3693

0.3659

0.3250

Trapezoidal

2.0283

1.1119

2.4544

Euler extrap.

2.3114

2.0172

1.9719

Semi-Trapezoidal extrap.

1.8687

1.9999

0.9834

Table 2.

Measured convergence rates Model 2.

Table 2.

Measured convergence rates Model 2.

Method

Call

Put

Indicator

Euler

0.4335

1.2898

0.8565

Semi-Trapezoidal

1.3025

0.7810

0.9518

FTE

1.2050

1.1733

1.0546

Symmetrized Euler

0.3028

0.3021

0.2421

Trapezoidal

1.8925

2.1272

1.6324

Euler extrap.

0.9483

1.4393

1.5966

Semi-Trapezoidal extrap.

1.4840

1.0481

1.2744

Table 3.

Measured convergence rates Model 3.

Table 3.

Measured convergence rates Model 3.

Method

Call

Put

Indicator

Euler

0.6977

0.5378

1.0695

Semi-Trapezoidal

1.6989

1.6551

1.6396

FTE

2.0091

1.7303

1.6008

Symmetrized Euler

1.0386

1.0426

0.9018

Trapezoidal

1.8682

1.6799

1.5219

Euler extrap.

1.1612

1.1857

2.2612

Semi-Trapezoidal extrap.

1.5660

1.0441

1.5979

Table 4.

Computational times (sec.) of the semi-exact schemes for 26 time steps and 2 × 107 paths.

Table 4.

Computational times (sec.) of the semi-exact schemes for 26 time steps and 2 × 107 paths.

Model

1

2

3

Euler

345.73

755.19

145.40

Semi-Trapezoidal

344.53

757.93

144.79

Trapezoidal

342.51

766.01

143.39

Euler extrap.

690.36

2335.94

307.62

Semi-Trapezoidal extrap.

686.55

2371.67

310.29

Table 5.

Computational times (sec.) of Euler-type discretizations for 26 time steps and 2 × 107 paths.

Table 5.

Computational times (sec.) of Euler-type discretizations for 26 time steps and 2 × 107 paths.

Model

1

2

3

FTE

142.3

138.37

141.53

Symmetrized Euler

141.64

140.98

141.67

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely

those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or

the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas,

methods, instructions or products referred to in the content.