The Weak Convergence Rate of Two Semi-Exact Discretization Schemes for the Heston Model

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 Eulerand 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.


Introduction and Main Results
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 S 0 , V 0 , κ, θ, σ > 0, µ ∈ R, ρ ∈ [−1, 1], T > 0 and independent Brownian motions W = (W t ) t∈[0,T] , B = (B t ) t∈ [0,T] , which are defined on a filtered probability space (Ω, F , (F t ) t∈[0,T] , P), 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 ν = 2κθ σ 2 will be an important parameter for our results. Throughout this article, the initial values S 0 , V 0 are assumed to be deterministic.
To price options with maturity at time T, one is interested in the value of where g : [0, ∞) → R is the payoff function. Closed formulae for E[g(S T )] are rarely known and often Monte Carlo methods are applied, for which in turn the simulation of S T 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 f : R → R with f (x) = g(exp(x)). While exact simulation schemes and their refinements are known (see, e.g., Broadie and Kaya (2006); Glasserman and Kim (2011); Malham and Wiese (2013); Smith (2007)), discretization schemes as, e.g., Altmayer and Neuenkirch (2017); Andersen (2008); Kahl and Jäckel (2006); Lord et al. (2009), are very popular for the Heston model. The latter discretization schemes can be easily extended to the multi-dimensional case and avoid computational bottlenecks of the exact schemes. In particular, Euler-type methods, such as the fully truncated Euler scheme, seem to be very efficient (see, e.g., Coskun and Korn ( 2018); Lord et al. (2009)), but no weak error analysis is available for them, up to the best of our knowledge.
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 V = (V t ) t∈[0,T] follows a non-central chi-square distribution, it can be simulated exactly. Trapezoidal discretizations of the first component X = (X t ) t∈ [0,T] lead to the trapezoidal scheme where 0 = t 0 < . . . < t k < . . . < t N = T, v k = V t k and ∆ k B = B t k +1 − B t k . 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 ν = 2κθ/σ 2 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 pol 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 .
In particular, for an equidistant discretization with t k = kT/N, k = 0, . . . , N, we have Here, u denotes the solution of the associated Kolmogorov PDE; see Equation (7).
In particular, for an equidistant discretization t k = kT/N, k = 0, . . . , N, 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
Remark 1. 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 E[ f (λ, X T , V T )] with λ ∈ [0, T], 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.
Remark 2. Property (6) allows to construct a second order scheme via extrapolation: If (6) holds, then Remark 3. 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.
Remark 4. 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 V = (V t ) t∈[0,T] 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 Z = (log(S t )) t∈[0,T] is discretized instead of the one for X = (X t ) t∈[0,T] .
Remark 5. 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.

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: The Feller index is ν = 2κθ σ 2 ≈ 0.63 in Model 1, ν ≈ 0.36 in Model 2, and ν ≈ 2.01 in Model 3. For each model, we use the following payoff functions:
Indicator: g 3 (S T ) = e −rT 1 [0,K] (S T ). 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 M = 2 · 10 7 independent copies g i (s (j) for each combination of model parameters, functional and number of steps N ∈ {2 1 , ..., 2 6 } where ∆t = T N . 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 p ref , i.e., e(N) = |p ref − p M,N |, and the error e(N) is plotted in Figures 1-18. We then measured the rate of convergence, i.e., the decay rate of e(N), 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 callprice formula given in Heston (1993) via the put-call-parity. The price of the digital option can be computed from the probability P 2 given in Heston (1993); it equals e rT (1 − P 2 ). 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).

Model 1
In 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.

Model 2
Here, we have an even lower Feller index of ν ≈ 0.36. 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 Figures 7,9,and 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.

Model 3
Here, we have the highest Feller Index with ν ≈ 2.01. 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 Figures 13, 15, and 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.

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 Tables 4 and 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.) Table 4. Computational times (sec.) of the semi-exact schemes for 2 6 time steps and 2 × 10 7 paths.

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.

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

Kolmogorov PDE
Recall that the stochastic integral equations for the log-Heston model for 0 ≤ s < t ≤ T 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 log(S t ) − log(S s ), since we have given below. To get the Kolmogorov backward PDE, we look at the following integral equations: We set and obtain for f : R × [0, ∞) → R 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: Moreover, we denote by |y| the standard Euclidean norm in R d . Let D ⊂ R d be a domain and q ∈ N. The set C q (D; R) is the set of all real-valued functions on D which are q-times continuously differentiable. For ε ∈ (0, 1), we denote by C q+ε (D; R) the set of all functions from C q (D; R) in which partial derivatives of order q are Höldercontinuous of order ε, and C q+ε c (D; R) is the set of all functions from C q+ε (D; R), who have compact support. Moreover, C q pol (D; R) is the set of functions g ∈ C q (D; R) such that there exist C, a > 0 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: Theorem 3. Let ε > 0 and f ∈ C 2+ε c (R × R + ; R). Then, there exists a constant c > 0, depending only on f , T, ρ, κ, θ 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 v ∈ [0, 1].
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 f x and dominated convergence we have whileũ fulfills the same PDE with terminal conditioñ Applying Theorem 3 now toû andũ, we obtain the following additional bounds (case (ii)) for the derivatives of u: c (R × R + ; R). Then, there exists a constant c > 0, depending only on f , T, ρ, κ, θ and σ such that the solution u of PDE (7) satisfies (ii) Let ε > 0 and f ∈ C 4+ε c (R × R + ; R). Then, there exists a constant c > 0, depending only on f , T, ρ, κ, θ 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: Theorem 4. Let q ∈ N, q ≥ 2 and suppose that f ∈ C 2q pol (R × R + ; R). Then, the solution u of PDE (7) satisfies u ∈ C q pol,T (R × R + ; R).
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.

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).
(2) For all p ≥ 1, there exist constants c > 0, depending only on p, κ, θ, σ, T, and V 0 , such that We will need the following bound on the growth of the L q -norm of a specific stochastic integral of the CIR process: Proof. With the Burkholder-Davis-Gundy inequality and the Hölder inequality, we have The assertion now follows from Lemma 1 (1).

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 (W t , B t ) t∈[0,T] a two-dimensional Brownian motion, then the Malliavin derivative measures the dependence of Y on (W, B). The Malliavin derivative is defined by a standard extension procedure: Let S be the set of smooth random variables of the form with ϕ ∈ C ∞ (R k ; R) bounded with bounded derivatives, h i ∈ L 2 ([0, T]; R 2 ), i = 1, . . . , k, and the stochastic integrals The derivative operator D of such a smooth random variable is defined as This operator is closable from L p (Ω) into L p Ω; H with H = L 2 ([0, T]; R 2 ) and the Sobolev space D 1,p denotes the closure of S with respect to the norm In particular, if D W 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 Y ∈ D 1,2 is independent of B, then D B Y = 0. For the CIR process, we will, therefore, have that D B V t = 0 for all t ∈ [0, T]. The derivative operator follows rules similar to ordinary calculus.
Proposition 1. Let X = (X 1 , ..., X d ) be a random variable with components in D 1,p . If then the chain rule holds: φ(X) ∈ D 1,p and For example, for a random variable Y ∈ D 1,p and g ∈ C 1 (R; R) 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 u ∈ L 2 Ω; L 2 ([0, T]; R 2 ) belongs to dom(δ), the domain of the divergence operator, then δ(u) is defined by the duality-also called integration by parts-relationship If u is adapted to the canonical filtration generated by (W, B) and satisfies E T 0 |u t | 2 dt < ∞, then u ∈ dom(δ) and δ(u) coincides with the Itō integral T 0 u 1 (s)dW s + T 0 u 2 (s)dB s . 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).

Lemma 3. Let t ∈ [0, T] and 2κθ
σ 2 > 1. Then, we have In Altmayer and Neuenkirch (2015), this and the integration by parts formula was used to establish under the assumption 2κθ σ 2 > 1 with G : R → R differentiable and g = G bounded, see Proposition 4.1 in Altmayer and Neuenkirch (2015). Indeed, t] (r) and the chain rule, i.e., 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 t s 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 (V, W) and B. Thus, we can conclude that

Properties of the Euler Discretization
Recall that the Euler discretization of the price process is given by We extend this discretization in every interval [t n , t n+1 ] as the following Itō process: Here, we have set n(t) := max{n ∈ {0, ..., N} : t n ≤ t}, η(t) := t n(t) and v k = V t k .
We have the following result on the Malliavin regularity of the Euler discretization: Lemma 4. Let t ∈ [0, T] and 2κθ σ 2 > 1. Then,x t ∈ D 1,∞ , and we have Proof. We havê Following the steps of the proof of Lemma 3.5 from Altmayer and Neuenkirch Altmayer and Neuenkirch (2017), we then havex t ∈ D 1,∞ exploiting that √ V t ∈ D 1,∞ and V t ∈ D 1,∞ under the assumption 2κθ σ 2 > 1. The chain rule from Proposition 1 yields Note that we write, in the following, v t instead of V t to unify the notation. With the above result, we can express E t η(t) √ v s dW s u xx (t,x t , v t ) without the second order derivative of u, which will be needed later on.

Lemma 5. Let t ∈ [0, T]. Under the assumptions of Theorem 3 and 2κθ
Proof. To avoid stronger restrictions on the Feller index we will use a localization procedure. So, for ε > 0, let ψ ε be a function such that 1.
Since (V, W) and B are independent, the chain rule from Proposition 1 implies Recall the integration by parts formula from Equation (8) where we now choose Before we can apply the integration by parts rule, we need to check whether for t > 0. We deduced these terms by using again the chain rule for D r Y. Note that the properties of the localizing function and Theorem 3 imply that are all uniformly bounded in (t, x, v). So, Equation (10) holds, then, due to Lemma 1, Lemma 3, Equation (9), and Lemma 4. Since t 0 1 √ v η(r) dB r is also well-defined by Lemma 1 due to 2κθ σ 2 > 1, we obtain now Due to Corollary 1 (i), not only u x but also u xx is bounded. Since ψ ε (v t ) → 1 almost surely for ε → 0 and |ψ ε (v t )| ≤ 1 for all ε > 0, the assertion follows now by dominated convergence using the Itô-isometry and again Lemma 1.
We also need the following L p -convergence result: Lemma 6. Let p ≥ 1. There exists a constant c > 0, depending only on p, T, ρ, κ, θ, σ and v 0 , such that sup

Proof.
We have Assume without loss of generality that p ≥ 2. Jensen's inequality and the Burkholder-Davis-Gundy inequality now imply that there exists a constant c > 0, depending only on p, T, the parameters of the CIR process, and v 0 , such that Since | √ x − √ y| ≤ |x − y| for x, y ≥ 0, the assertion follows from Lemma 1.
Straightforward calculations also yield the following L p -smoothness result for the Euler-type scheme: Lemma 7. Let p ≥ 1. There exists a constant c > 0, depending only on p, T, ρ, κ, θ, σ, and v 0 , such that E|x t −x s | p ≤ c · |t − s| p/2 for all s, t ∈ [0, T].

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 Here, we have set again n(t) := max{n ∈ {0, ..., N} : t n ≤ t}, η(t) := t n(t) , v k = V t k , and we also write again v t , instead of V t , to unify the notation.
We need the following result on the Malliavin regularity of the semi-trapezoidal scheme: Lemma 8. Let t ∈ [0, T] and 2κθ σ 2 > 1. Then, we havex t ∈ D 1,∞ and Proof. We already know that v t ∈ D 1,∞ and √ v t ∈ D 1,∞ . We can writex t aŝ Following the steps of the proof of Lemma 3.5 from Altmayer and Neuenkirch (2017), we then also havex t ∈ D 1,∞ . The chain rule from Proposition 1 yields 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: Lemma 9. Let t ∈ [0, T]. Under the assumptions of Theorem 3 and 2κθ By similar calculations as for the Euler scheme, we also have: Lemma 10. Let p ≥ 1. There exists a constant c > 0, depending only on p, T, ρ, κ, θ, σ, and v 0 , such that sup Lemma 11. Let p ≥ 1. There exists a constant c > 0, depending only on p, T, ρ, κ, θ, σ, and v 0 , such that E|x t −x s | p ≤ c · |t − s| p/2 for all s, t ∈ [0, T].

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 f , T, ρ, κ, θ, σ, and v 0 , x 0 , will be denoted by c, regardless of their value.

The Euler Scheme: Expanding the Error
, the weak error is a telescoping sum of local errors: With the Itō formula and the Kolmogorov backward PDE evaluated at (t,x t , v t ), we obtain we have e n = e (1) n with e (1) By Theorem 3 and Corollary 1, we have that u x and u xx are bounded. So, Lemma 1 implies that Moreover, with the law of total expectation and the adaptedness ofx η(t) and v η(t) , we have due to the martingale property of the Itō integral. Therefore, we can write and obtain e (1) In the same way, we have e (2) Summarizing this preliminary part, we have obtained wherẽ e (1) 4.2. The Euler Scheme: Case (i)

So, it remains to analyzeẽ
(1) n andẽ (2) n under the regularity of Theorem 3 (i). We start withẽ (1) n . The mean value theorem and where we used With the Minkowski inequality and Lemma 1, it holds that Finally, with the Minkowski inequality, the Burkholder-Davis-Gundy inequality, Lemma 1, and Lemma 7, we obtain for all p ≥ 1 that where c is in particular independent of t ∈ [0, T]. The Hölder inequality then gives e (1) n = O((∆t) 2 ).

Forẽ
(2) n , 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 1/2. First, recall that, by Lemma 5, we have Moreover, note that we also have and recall that Thus, we can writẽ e (2) with Before we analyze this expression further, it remains to show (14). Using the law of total expectation, the adaptedness ofx η(t) , v η(t) and of the Itō integrals, we have due to the properties of the Itō integral, Equation (14) follows.
Using the mean value theorem in Equation (15), we obtain Therefore, By Lemmas 1 and 2, it holds that σ 2 ) and p ∈ [0, 2κθ σ 2 ). So, the Hölder inequality leads to 3σ 2 − 1) and c in particular independent of t ∈ [0, T]. With the Cauchy-Schwarz, Burkholder-Davis-Gundy, and Minkowski inequalities for p ≥ 1, it follows that With the Hölder inequality, we now have Therefore, which concludes the proof of this part.

The Euler Scheme: Case (ii)
Starting from Equation (11) and using now the bounds of Corollary 1 for u xx , u xv , u xxx , and u xxv , the assertion follows from a direct application of the mean value theorem to (12) and (13), together with the Lemmata 1 and 7.

Semi-Trapezoidal Rule: Expanding the Error
We look again at the telescoping sum of local errors With the Itō formula and the Kolmogorov backward PDE evaluated at (t,x t , v t ), we have we obtain e n = e n + e n + e (3) n with e (1) using Lemma 1. Moreover, since u x and u xx are bounded by Theorem 3 and Corollary 1 (i), we obtain similar to the calculations for the Euler scheme that e (1) withẽ (1) √ v s dW s u xx (t,x t , v t ) − u xx (t,x η(t) , v η(t) ) dt.
4.5. Semi-Trapezoidal Rule: Case (i) Since Lemma 8 givesx t ∈ D 1,∞ and we can proceed here in the same way as for the Euler scheme by using the Lemmata 9 and 11.

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.

Proof of Theorem 2
Now, we derive the error expansion under the regularity of Theorem 4 with q = 4, i.e., we have u ∈ C 4 pol,T (R × R + ; R). E|X t | p < ∞ and E|X t − X s | p ≤ c · |t − s| p/2 , s, t ∈ [0, T], for all p ≥ 1. We will use this in the following at several places without explicitly mentioning it. Recall that we obtained e (1)
If higher derivatives of u are available, then we can analyze E t η(t) √ v s dW s u x (t,x t , v t ) and E t η(t) √ v s dW s u xx (t,x t , v t ) via another application of Itō's lemma. So, let k : [0, T] × R × [0, ∞) → R be a C 1,2 -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 k x and k v have polynomial growth, then an application of the Itō isometry and the martingale property of the Itō integral yield E|X t | p < ∞ and E|X t − X s | p ≤ c · |t − s| p/2 , s, t ∈ [0, T], for all p ≥ 1. We will use this in the following at several places without explicitly mentioning it.

Semi-Trapezoidal Scheme: Preliminaries
We will now take also a closer look at the error of the semi-trapezoidal discretization for u ∈ C 4 pol,T (R × R + ; R). Recall that and so where H(s, t,x s ,x t , v s , v t ) = − (1 − ρ 2 ) 2 κ(θ − v s )u xx (t,x t , v t ) + σ 2 v s u xxv (s,x s , v s ) .
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 s, t ∈ [t n , t n+1 ] that E[H(s, t,x s ,x t , v s , v t )] = E[H(t n , t n , X t n , X t n , V t n , V t n )] + O((∆t) 1/4 ).
In particular, for an equidistant discretization t k = kT/N, k = 0, . . . , N, we have E H(t n , t n , X t n , X t n , V t n , V t n )∆t + O((∆t) 5/4 ) and the convergence N−1 ∑ n=0 E H(t n , t n , X t n , X t n , V t n , V t n )∆t → T 0 E H(t, t, X t , X t , V t , V t )dt for ∆t → 0 concludes the proof of Theorem 2 (ii).