The weak convergence rate of two semi-exact discretization schemes for the Heston model
Abstract
EconStor is a publication server for scholarly economic literature, provided as a non-commercial public service by the ZBW.
Full text
Mickel, Annalena; Neuenkirch, Andreas Article The weak convergence rate of two semi-exact discretization schemes for the Heston model Risks Provided in Cooperation with: MDPI – Multidisciplinary Digital Publishing Institute, Basel Suggested Citation: Mickel, Annalena; Neuenkirch, Andreas (2021) : The weak convergence rate of two semi-exact discretization schemes for the Heston model, Risks, ISSN 2227-9091, MDPI, Basel, Vol. 9, Iss. 1, pp. 1-38, https://doi.org/10.3390/risks9010023 This Version is available at: https://hdl.handle.net/10419/258113 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. https://creativecommons.org/licenses/by/4.0/
risks Article The Weak Convergence Rate of Two Semi-Exact Discretization Schemes for the Heston Model Annalena Mickel 1,2 and Andreas Neuenkirch 2,* Citation: Mickel, Annalena, and Andreas Neuenkirc. 2021. The Weak Convergence Rate of two Semi-Exact Discretization Schemes for the Heston Model. Risks 9: 23. https://doi.org/ 10.3390/risks9010023 Received: 15 October 2020 Accepted: 28 December 2020 Published: 12 January 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: c 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). 1DFG Research Training Group 1953, University of Mannheim, B6, 26, D-68131 Mannheim, Germany; [email protected] 2Mathematical Institute, University of Mannheim, B6, 26, D-68131 Mannheim, Germany *Correspondence: [email protected] Abstract: 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. Keywords: Heston model; discretization schemes for SDEs; exact simulation of the CIR process; Kolmogorov PDE; Malliavin calculus MSC: 60H07; 60H35; 65C05; 91G60 1. 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 Sand its volatility V: dSt=µStdt +√VtStρdWt+q1−ρ2dBt, dVt=κ(θ−Vt)dt +σ√VtdWt, (1) with S0 , V0 , κ , θ , σ> 0, µ∈R , ρ∈[− 1,1 ] , T> 0 and independent Brownian motions W= (Wt)t∈[0,T] , B= (Bt)t∈[0,T] , which are defined on a filtered probability space (Ω , F , (Ft)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 S0,V0are assumed to be deterministic. To price options with maturity at time T, one is interested in the value of E[g(ST)], where g:[ 0, ∞)→R is the payoff function. Closed formulae for E[g(ST)] are rarely known and often Monte Carlo methods are applied, for which in turn the simulation of ST Risks 2021,9, 23. https://doi.org/10.3390/risks9010023 https://www.mdpi.com/journal/risks
Risks 2021,9, 23 2 of 38 is required. Usually, the log-Heston model instead of the Heston model is considered in numerical practice. This yields the SDE d(log(St)) = µ−1 2Vtdt +√VtdρWt+q1−ρ2Bt, dVt=κ(θ−Vt)dt +σ√VtdWt, (2) 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 dXt=ρκ σ−1 2Vtdt +q1−ρ2√VtdBt, dVt=κ(θ−Vt)dt +σ√VtdWt. (3) Moreover, since the transition density of the CIR process V= (Vt)t∈[0,T] follows a non-central chi-square distribution, it can be simulated exactly. Trapezoidal discretizations of the first component X= (Xt)t∈[0,T]lead to the trapezoidal scheme xk+1=xk+ρκ σ−1 2vk+1+vk 2(tk+1−tk) +q1−ρ2rvk+1+vk 2∆kB,k=0, . . . , N−1, where 0 =t0<. . . <tk<. . . <tN=T , vk=Vtk and ∆kB=Btk+1−Btk . 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 ¯ o-Taylor expansions turned out to be not feasible. So, instead, we analyzed the following two semi-exact discretization schemes: the Euler-type scheme xk+1=xk+ρκ σ−1 2vk(tk+1−tk) + q1−ρ2√vk∆kB(4) and the semi-trapezoidal scheme xk+1=xk+ρκ σ−1 2vk+1+vk 2(tk+1−tk) + q1−ρ2√vk∆kB. (5) 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
Risks 2021,9, 23 3 of 38 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 ∆t=max k=1,...,N|tk−tk−1| 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. Theorem 1. Let ε>0. (i) If f ∈C2+ε c(R×R+;R)and 2κθ σ2>3 2, then both schemes satisfy E[f(xN,vN)]−E[f(XT,VT)]=O(∆t). (ii) If f ∈C4+ε c(R×R+;R), then both schemes satisfy E[f(xN,vN)]−E[f(XT,VT)]=O(∆t). Assuming more smoothness of f, we obtain more detailed results: Theorem 2. Suppose that f ∈C8 pol(R×R+;R). (i) Then, the Euler scheme (4)satisfies E[f(xN,vN)]−E[f(XT,VT)]= N−1 ∑ n=0Ztn+1 tnZt tn E[H(s,t,ˆ xs,ˆ xt,Vs,Vt)]dsdt +O((∆t)2), where H(s,t,ˆ xs,ˆ xt,Vs,Vt) = 1 2−ρκ σκ(θ−Vs)ux(t,ˆ xt,Vt) + σ2Vsuxv(s,ˆ xs,Vs) −(1−ρ2) 2κ(θ−Vs)uxx(t,ˆ xt,Vt) + σ2Vsuxxv(s,ˆ xs,Vs) and ˆ xt=xn+ρκ σ−1 2vn(t−tn) + q1−ρ2√vn(Bt−Btn),t∈[tn,tn+1], for n =0, . . . , N−1. In particular, for an equidistant discretization with tk=kT/N, k =0, . . . , N, we have lim N→∞N(E[f(xN,vN)]−E[f(XT,VT)]) =T 2ZT 0 E[H(t,t,Xt,Xt,Vt,Vt)]dt. Here, u denotes the solution of the associated Kolmogorov PDE; see Equation (7).
Risks 2021,9, 23 4 of 38 (ii) For the semi-trapezoidal scheme (5), we have E[f(xN,vN)]−E[f(XT,VT)]= N−1 ∑ n=0Ztn+1 tnZt tn E[H(s,t,ˆ xs,ˆ xt,Vs,Vt)]dsdt +O((∆t)2), where H(s,t,ˆ xs,ˆ xt,Vs,Vt) = −(1−ρ2) 2κ(θ−Vs)uxx(t,ˆ xt,Vt) + σ2Vsuxxv(s,ˆ xs,Vs) and ˆ xt=xn+ρκ σ−1 2Vt+vn 2(t−tn) + q1−ρ2√vn(Bt−Btn),t∈[tn,tn+1], for n =0, . . . , N−1. In particular, for an equidistant discretization tk=kT/N, k =0, . . . , N, it holds lim N→∞N(E[f(xN,vN)]−E[f(XT,VT)]) =T 2ZT 0 E[H(t,t,Xt,Xt,Vt,Vt)]dt. 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 E[f(xN,vN)]−E[f(XT,VT)]=T 2ZT 0 E[H(t,t,xt,xt,vt,vt)]dt·N−1+O(N−2). (6) However, to establish this, we would require error estimates for functionals of the type E[f(λ , XT , VT)] 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 YN=2f(x2N,v2N)−f(xN,vN) satisfies EYN=Ef(XT,VT) + O((∆t)2), where (x2N,v2N)uses the stepsize T/(2N)and (xN,vN)the stepsize T/N. 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
Risks 2021,9, 23 5 of 38 studied in Altmayer and Neuenkirch (2017); Lord et al. (2009). For example, the so-called absorbed Euler discretization reads as zk+1=zk−1 2vk(tk+1−tk) + √vkρWtk+1−Wtk+q1−ρ2Btk+1−Btk, vk+1=vk+κ(θ−vk)(tk+1−tk) + σ√vkWtk+1−Wtk+. Here, the volatility V= (Vt)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(St))t∈[0,T]is discretized instead of the one for X = (Xt)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. 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 1and 2. We use the following model parameters: Model 1: S0= 100, V0= 0.010201, K= 100, κ= 6.21, θ= 0.019, σ= 0.61, ρ=− 0.7, T=1,r=0.0319; Model 2: S0= 100, V0= 0.09, K= 100, κ= 2, θ= 0.09, σ= 1, ρ=− 0.3, T= 5, r= 0.05; Model 3: S0= 100, V0= 0.0457, K= 100, κ= 5.07, θ= 0.0457, σ= 0.48, ρ=− 0.767, T=2,r=0.00. 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: 1. European Call: g1(ST) = e−rT max{ST−K,0}; 2. European Put: g2(ST) = e−rT max{K−ST,0}; 3. Indicator: g3(ST) = e−rT1[0,K](ST). 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 gi(s(j) N),j=1, . . . , M, of gi(sN)to estimate E(gi(sN)) by pM,N=1 M M ∑ j=1 gi(s(j) N) 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 pref , i.e., e(N) = |pref −pM,N|,
Risks 2021,9, 23 6 of 38 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 P2 given in Heston (1993); it equals erT(1−P2) . 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). 2.1. Model 1 In Table 1, we can see the measured convergence rates for this model with a Feller index of ν≈0.63. The associated plots are shown in Figures 1–6. 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 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 1and 2could 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.
Risks 2021,9, 23 7 of 38 Figure 1. Call Model 1. Figure 2. Call Model 1. Figure 3. Put Model 1.
Risks 2021,9, 23 8 of 38 Figure 4. Put Model 1. Figure 5. Indicator Model 1. Figure 6. Indicator Model 1. 2.2. 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
Risks 2021,9, 23 15 of 38 ut(t,x,v) = −ρκ σ−1 2vux(t,x,v)−κ(θ−v)uv(t,x,v) −v 21−ρ2uxx(t,x,v) + σ2uvv(t,x,v),t∈(0, T),x∈R,v>0, u(T,x,v) = f(x,v),x∈R,v≥0. (7) 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 l= (l1 ,..., ld)∈Nd , we define |l|=∑d j=1lj and for y∈Rd , we define ∂l y=∂l1 y1···∂ld yd . Moreover, we denote by |y| the standard Euclidean norm in Rd . Let D ⊂ Rd be a domain and q∈N . The set Cq(D;R) is the set of all real-valued functions on D which are q -times continuously differentiable. For ε∈( 0,1 ) , we denote by Cq+ε(D;R) the set of all functions from Cq(D;R) in which partial derivatives of order q are Höldercontinuous of order ε , and Cq+ε c(D;R) is the set of all functions from Cq+ε(D;R) , who have compact support. Moreover, Cq pol(D;R) is the set of functions g∈Cq(D;R) such that there exist C,a>0 for which |∂l yg(y)| ≤ C(1+|y|a),y∈ D,|l| ≤ q. Finally, we denote by Cq pol,T(D;R) the set of functions v∈Cbq/2c,q pol ([0, T)×D;R) such that there exist C,a>0 for which sup t<T|∂k t∂l yv(t,y)| ≤ C(1+|y|a),y∈ D,2k+|l| ≤ q. 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 ε> 0and f∈C2+ε 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 sup (t,x,v)∈[0,T]×R×[0,∞)|u(t,x,v)|+|∂tu(t,x,v)|+|∂vu(t,x,v)|+|∂xu(t,x,v)|≤c, sup (t,x,v)∈[0,T]×R×[0,1]|v∂xxu(t,x,v)|+|v∂xvu(t,x,v)|+|v∂vvu(t,x,v)|≤c, sup (t,x,v)∈[0,T]×R×[1,∞)|∂xxu(t,x,v)|+|∂xvu(t,x,v)|+|∂vvu(t,x,v)|≤c. 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 vfor 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 ˆ u(t,x,v):=ux(t,x,v) = Ehfx(Xt,x,v T,Vt,v T)i, ˜ u(t,x,v):=uxx(t,x,v) = Ehfxx(Xt,x,v T,Vt,v T)i.
Risks 2021,9, 23 16 of 38 This is well defined: by continuity and boundedness of fx and dominated convergence we have ux(t,x,v) = ∂ ∂xEhf(Xt,x,v T,Vt,v T)i=∂ ∂xEhf(x+Zt,v T,Vt,v T)i =lim δ→0 E"f(x+δ+Zt,v T,Vt,v T)−f(x+Zt,v T,Vt,v T) δ# =lim δ→0Z1 0 Ehfx(x+δλ +Zt,v T,Vt,v T)idλ =Z1 0 Ehfx(x+Zt,v T,Vt,v T)idλ =Ehfx(x+Zt,v T,Vt,v T)i=Ehfx(Xt,x,v T,Vt,v T)i with Zt,v T=ρκ σ−1 2Zt sVs,v rdr +q1−ρ2Zt sqVs,v rdBr. An analogous calculation for uxx(t , x , v) shows that uxx(t , x , v) = Ehfxx(Xt,x,v T,Vt,v T)i . Thus, uxx is also bounded, if f∈C2+ε c(R×R+ ; R) . Moreover, ˆ u fulfills the Kolmogorov backward PDE ˆ ut(t,x,v) = −ρκ σ−1 2vˆ ux(t,x,v)−κ(θ−v)ˆ uv(t,x,v) −v 21−ρ2ˆ uxx(t,x,v) + σ2ˆ uvv(t,x,v),t∈(0, T),x∈R,v>0, ˆ u(T,x,v) = fx(x,v),x∈R,v≥0, while ˜ ufulfills the same PDE with terminal condition ˜ u(T,x,v) = fxx(x,v),x∈R,v≥0. Applying Theorem 3now to ˆ u and ˜ u , we obtain the following additional bounds (case (ii)) for the derivatives of u: Corollary 1. (i) Let ε> 0and f∈C2+ε 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 sup (t,x,v)∈[0,T]×R×[0,∞)|∂xxu(t,x,v)| ≤ c. (ii) Let ε> 0and f∈C4+ε 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 sup (t,x,v)∈[0,T]×R×[0,∞)|∂xvu(t,x,v)|+|∂xxu(t,x,v)|+|∂xxvu(t,x,v)|+|∂xxxu(t,x,v)|≤c. 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≥ 2and suppose that f∈C2q pol(R×R+ ; R) . Then, the solution u of PDE (7)satisfies u ∈Cq pol,T(R×R+;R).
Risks 2021,9, 23 17 of 38 In contrast to the results of Feehan and Pop, the result of Briani et al. requires more smoothness of fbut 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). Lemma 1. (1) We have E"sup t∈[0,T] Vp t#<∞ for all p ≥1and sup t∈[0,T] EVp t<∞iff p >−2κθ σ2. (2) For all p ≥1, there exist constants c >0, depending only on p,κ,θ,σ,T,and V0, such that E|Vt−Vs|p≤c·|t−s|p/2,s,t∈[0, T]. We will need the following bound on the growth of the Lq -norm of a specific stochastic integral of the CIR process: Lemma 2. For all q ∈h2, 4κθ σ2, it holds that sup t∈[0,T] t−q/2 EZt 0 1 √Vu dBu q<∞. Proof. With the Burkholder-Davis-Gundy inequality and the Hölder inequality, we have t−q/2EZt 0 1 √Vu dBu q≤t−q/2E"Zt 0 1 Vu du q/2# ≤t−q/2E Zt 01 Vuq/2 du 2/qZt 0dr (q−2)/q!q/2 =t−q/2Zt 0 E"1 Vuq/2#dut(q−2)/2 ≤sup u∈[0,t] EhV−q/2 ui for all t∈[0, T]. 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 (Wt , Bt)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 S=ϕZT 0h1(s)d(Ws,Bs), . . . , ZT 0hk(s)d(Ws,Bs)
Risks 2021,9, 23 18 of 38 with ϕ∈C∞(Rk ; R) bounded with bounded derivatives, hi∈L2([ 0, T] ; R2) , i= 1, . . . , k , and the stochastic integrals ZT 0hj(s)d(Ws,Bs) = ZT 0h(1) j(s)dWs+ZT 0h(2) j(s)dBs. The derivative operator Dof such a smooth random variable is defined as DS = k ∑ i=1 ∂ϕ ∂xiZT 0h1(s)d(Ws,Bs), . . . , ZT 0hk(s)d(Ws,Bs)hi. This operator is closable from Lp(Ω) into LpΩ ; H with H=L2([ 0, T] ; R2) and the Sobolev space D1,pdenotes the closure of Swith respect to the norm kYk1,p=E|Y|p+EZT 0|DsY|2ds p1/p . In particular, if DW denotes the first component of the Malliavin derivative, i.e., the derivative with respect to W, we have DW tY=1[0,t]if Y=W 0 if Y=B and vice versa for the derivative with respect to B, i.e., DB tY=1[0,t]if Y=B 0 if Y=W This, in particular, implies that, if Y∈D1,2 is independent of B, then DBY=0. For the CIR process, we will, therefore, have that DBVt=0 for all t∈[0, T]. The derivative operator follows rules similar to ordinary calculus. Proposition 1. Let X = (X1,..., Xd)be a random variable with components in D1,p. If (i) φ:Rd→Ris in C1(Rd;R), (ii) φ(X)∈Lp(Ω), (iii) ∂iφ(X)·DXi∈Lp(Ω;H)for all i =1, ..., d, then the chain rule holds: φ(X)∈D1,pand Dφ(X) = d ∑ i=1 ∂iφ(X)·DXi. For example, for a random variable Y∈D1,p and g∈C1(R ; R) with bounded derivative, the chain rule reads as Dg(Y) = g0(Y)DY. Another simple example for the application of this chain rule is DW r(Wt−Ws)2=2(Wt−Ws)1(s,t](r),r,s,t∈[0, T],s≤t. The divergence operator δ is the adjoint of the derivative operator. If a random variable u∈L2Ω ; L2([ 0, T] ; R2) belongs to dom(δ) , the domain of the divergence operator, then δ(u)is defined by the duality—also called integration by parts—relationship E[Yδ(u)] = EZT 0hDsY,usidsfor all Y∈D1,2. (8)
Risks 2021,9, 23 19 of 38 If u is adapted to the canonical filtration generated by (W , B) and satisfies ERT 0|ut|2dt <∞ , then u∈dom(δ)and δ(u)coincides with the It¯ o integral RT 0u1(s)dWs+RT 0u2(s)dBs. For the Malliavin regularity of the CIR process, the following is well known. See, e.g., Proposition 4.5 and Theorem4.6 in Altmayer (2015) orProposition4.1 in Alos and Ewald (2008) . Lemma 3. Let t ∈[0, T]and 2κθ σ2>1. Then, we have √Vt∈D1,∞and Vt∈D1,∞with Dr√Vt=σ 2expZt rσ2 8−κθ 21 Vu−κ 2du1[0,t](r),r∈[0, T]. In Altmayer and Neuenkirch (2015), this and the integration by parts formula was used to establish E(g(XT)) = 1 Tp1−ρ2·EG(XT)·ZT 0 1 √Vt dBt, under the assumption 2κθ σ2> 1 with G:R→R differentiable and g=G0 bounded, see Proposition 4.1 in Altmayer and Neuenkirch (2015). Indeed, using ut= 1 /√Vt and ERT 0|ut|2dt <∞,DB rXt=p1−ρ2√Vt1[0,t](r)and the chain rule, i.e., DB rG(XT) = g(XT)DB rXT, we have EG(XT)·ZT 0 1 √Vt dBt=EZT 0g(XT)·DB tXT·1 √Vt dt =EZT 0g(XT)q1−ρ2√Vt·1 √Vt dt=Tq1−ρ2E[g(XT)], where the first equality is due to the integration by parts formula. In Lemmas 5and 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 Rt s√VudWu. Since Zt s √VudWu=1 σVt−Vs−κθ(t−s) + κZt sVudu, we obtain DW rZt s √VudWu=1 σDW r(Vt−Vs) + κZt sDW rVudu, DB rZt s √VudWu=0, by exchanging the Riemann integral and the Malliavin derivative (via a standard approximation argument for the Riemann integral, Lemma 3and Lemma 1.2.3 in Nualart (1995)) and the independence of (V,W)and B. Thus, we can conclude that Zt s √VudWu∈D1,∞, 0 ≤s<t≤T. (9) 3.4. Properties of the Euler Discretization Recall that the Euler discretization of the price process is given by xk+1=xk+ρκ σ−1 2vk(tk+1−tk) + q1−ρ2√vk∆kB
Risks 2021,9, 23 20 of 38 with ∆kB=Btk+1−Btk . We extend this discretization in every interval [tn , tn+1] as the following It¯ o process: ˆ xt=xn(t)+ρκ σ−1 2Zt η(t)vn(t)ds +q1−ρ2Zt η(t)pvn(t)dBs. Here, we have set n(t):=max{n∈ {0,..., N}:tn≤t},η(t):=tn(t)and vk=Vtk. We have the following result on the Malliavin regularity of the Euler discretization: Lemma 4. Let t ∈[0, T]and 2κθ σ2>1. Then, ˆ xt∈D1,∞, and we have DB rˆ xt=q1−ρ2pvn(r)1[0,t](r). Proof. We have ˆ xt=ˆ xη(t)+ρκ σ−1 2vn(t)(t−η(t)) + q1−ρ2pvn(t)(Bt−Bη(t)) and ˆ xη(t)=ρκ σ−1 2n(t)−1 ∑ k=0 vk(tk+1−tk) + q1−ρ2 n(t)−1 ∑ k=0 √vk(Btk+1−Btk). Following the steps of the proof of Lemma 3.5 from Altmayer and Neuenkirch Altmayer and Neuenkirch (2017), we then have ˆ xt∈D1,∞ exploiting that √Vt∈D1,∞ and Vt∈D1,∞under the assumption 2κθ σ2>1. The chain rule from Proposition 1yields DB rˆ xη(t)=q1−ρ2 n(t)−1 ∑ k=0 √vk1(tk,tk+1](r) and DB rˆ xt=DB rˆ xη(t)+q1−ρ2pvn(t)1(tn(t),t](r). Note that we write, in the following, vt instead of Vt to unify the notation. With the above result, we can express EhRt η(t)√vsdWsuxx(t,ˆ xt,vt)i without the second order derivative of u, which will be needed later on. Lemma 5. Let t ∈[0, T]. Under the assumptions of Theorem 3and 2κθ σ2>1, we have EZt η(t)√vsdWsuxx(t,ˆ xt,vt)=1 tp1−ρ2E"Zt η(t)√vsdWsux(t,ˆ xt,vt)Zt 0 1 √vη(r) dBr#. 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. ψε:R→Ris continuously differentiable with bounded derivative, 2. 0 ≤ψε(x)≤1 on [0, ∞), 3. ψε(x) = 1 on [2ε,∞), 4. ψε(x) = 0 on (−∞,ε].
Risks 2021,9, 23 21 of 38 Since (V,W)and Bare independent, the chain rule from Proposition 1implies DB rZt η(t)√vsdWsψε(vt)ux(t,ˆ xt,vt)=Zt η(t)√vsdWsψε(vt)uxx(t,ˆ xt,vt)DB rˆ xt with DB rˆ xt=p1−ρ2√vη(r)1[0,t](r) . Recall the integration by parts formula from Equation (8), i.e., EYZT 0u1(s)dWs+ZT 0u2(s)dBs=EZT 0hDsY,usids, where we now choose DrY= DW rRt η(t)√vsdWsψε(vt)ux(t,ˆ xt,vt) DB rRt η(t)√vsdWsψε(vt)ux(t,ˆ xt,vt) ,ur= 0 1 √vη(r)1[0,t](r)!. Before we can apply the integration by parts rule, we need to check whether ZT 0 E" DW rZt η(t)√vsdWsψε(vt)ux(t,ˆ xt,vt) 2#dr <∞, ZT 0 E"Zt η(t)√vsdWsψ0 ε(vt)ux(t,ˆ xt,vt)DW rvt 2#dr <∞, ZT 0 E"Zt η(t)√vsdWsψε(vt)uxx(t,ˆ xt,vt)DW rˆ xt+uxv(t,ˆ xt,vt)DW rvt 2#dr <∞, ZT 0 E"Zt η(t)√vsdWsψε(vt)uxx(t,ˆ xt,vt)DB rˆ xt 2#dr <∞, (10) for t> 0. We deduced these terms by using again the chain rule for DrY . Note that the properties of the localizing function and Theorem 3imply that ψε(v)ux(t,x,v),ψ0 ε(v)ux(t,x,v),ψε(v)uxx(t,x,v),ψε(v)uxv(t,x,v) 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 Rt 01 √vη(r)dBris also well-defined by Lemma 1due to 2κθ σ2>1, we obtain now EZt η(t)√vsdWsψε(vt)uxx(t,ˆ xt,vt) =1 tEZt 0Zt η(t)√vsdWsψε(vt)uxx(t,ˆ xt,vt)dr =1 tp1−ρ2E"Zt 0DB rZt η(t)√vsdWsψε(vt)ux(t,ˆ xt,vt)1 √vη(r) dr# =1 tp1−ρ2E"Zt η(t)√vsdWsψε(vt)ux(t,ˆ xt,vt)Zt 0 1 √vη(r) dBr#.
Risks 2021,9, 23 22 of 38 Due to Corollary 1(i), not only ux but also uxx is bounded. Since ψε(vt)→ 1 almost surely for ε→ 0 and |ψε(vt)| ≤ 1 for all ε> 0, the assertion follows now by dominated convergence using the Itô-isometry and again Lemma 1. We also need the following Lp-convergence result: Lemma 6. Let p≥ 1. There exists a constant c> 0, depending only on p , T , ρ , κ , θ , σ and v0 , such that sup t∈[0,T] E|Xt−ˆ xt|p≤c·(∆t)p/4. Proof. We have Xt−ˆ xt=ρκ σ−1 2Zt 0(vu−vη(u))du +q1−ρ2Zt 0(√vu−pvη(u))dBu. 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 v0, such that E|Xt−ˆ xt|p≤cZt 0 E|vu−vη(u)|pdu +cZt 0 E|√vu−pvη(u)|pdu. Since |√x−√y| ≤ p|x−y|for x,y≥0, the assertion follows from Lemma 1. Straightforward calculations also yield the following Lp -smoothness result for the Euler-type scheme: Lemma 7. Let p≥ 1. There exists a constant c> 0, depending only on p , T , ρ , κ , θ , σ , and v0 , such that E|ˆ xt−ˆ xs|p≤c·|t−s|p/2 for all s,t∈[0, T]. 3.5. Properties of the Semi-Trapezoidal Rule Recall that our semi-trapezoidal rule reads as xk+1=xk+ρκ σ−1 21 2(vk+1+vk)(tk+1−tk) + q1−ρ2√vk∆kB =xk+ρκ σ−1 2vk(tk+1−tk) + q1−ρ2√vk∆kB +ρκ σ−1 21 2(vk+1−vk)(tk+1−tk). Again, we write the scheme as a time-continuous process: ˆ xt=xn(t)+ρκ σ−1 2vn(t)(t−η(t)) + q1−ρ2pvn(t)Bt−Bη(t) +ρκ σ−1 21 2vt−vn(t)(t−η(t)). Expanding the last term with It¯ o’s lemma, we obtain ˆ xt=xn(t)+Zt η(t)asds +Zt η(t)bsdBs+Zt η(t)csdWs
Risks 2021,9, 23 23 of 38 with at:=ρκ σ−1 2vn(t)+1 2(t−η(t))κ(θ−vt) + 1 2(vt−vn(t)), bt:=q1−ρ2pvn(t), ct:=ρκ σ−1 21 2(t−η(t))σ√vt. Here, we have set again n(t):=max{n∈ { 0,..., N}:tn≤t} , η(t):=tn(t) , vk=Vtk , and we also write again vt, instead of Vt, 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 have ˆ xt∈D1,∞and DB rˆ xt=q1−ρ2pvη(r)1[0,t](r). Proof. We already know that vt∈D1,∞and √vt∈D1,∞. We can write ˆ xtas ˆ xt=ˆ xη(t)+1 2ρκ σ−1 2vt+vη(t)(t−η(t)) + q1−ρ2Zt η(t)pvη(t)dBs with ˆ xη(t)=1 2ρκ σ−1 2n(t)−1 ∑ k=0 (vk+1+vk)(tk+1−tk) + q1−ρ2 n(t)−1 ∑ k=0 √vk(Btk+1−Btk). Following the steps of the proof of Lemma 3.5 from Altmayer and Neuenkirch (2017), we then also have ˆ xt∈D1,∞. The chain rule from Proposition 1yields DB rˆ xη(t)=q1−ρ2 n(t)−1 ∑ k=0 √vk1(tk,tk+1](r) and DB rˆ xt=DB rˆ xη(t)+q1−ρ2pvn(t)1(tn(t),t](r). 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 3and 2κθ σ2>1, we have EZt η(t)√vsdWsuxx(t,ˆ xt,vt)=1 tp1−ρ2E"Zt η(t)√vsdWsux(t,ˆ xt,vt)Zt 0 1 √vη(r) dBr#. 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 v0 , such that sup t∈[0,T] E|Xt−ˆ xt|p≤c·(∆t)p/4.
Risks 2021,9, 23 24 of 38 Lemma 11. Let p≥ 1. There exists a constant c> 0, depending only on p , T , ρ , κ , θ , σ , and v0 , such that E|ˆ xt−ˆ xs|p≤c·|t−s|p/2 for all s,t∈[0, T]. 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 ∆t=max k=1,...,N|tk−tk−1|, and depend only f,T,ρ,κ,θ,σ, and v0,x0, will be denoted by c, regardless of their value. 4.1. The Euler Scheme: Expanding the Error Since u(T , xN , vN) = Ef(xN , vN) and u( 0, x0 , v0) = Ef(XT , VT) , the weak error is a telescoping sum of local errors: |E[f(xN,vN)]−E[f(XT,VT)]|= N ∑ n=1 E[u(tn,xn,vn)−u(tn−1,xn−1,vn−1)] . With the It ¯ o formula and the Kolmogorov backward PDE evaluated at (t , ˆ xt , vt) , we obtain en:=E[u(tn+1,xn+1,vn+1)−u(tn,xn,vn)] =Ztn+1 tn Eut(t,ˆ xt,vt) + ρκ σ−1 2vn(t)ux(t,ˆ xt,vt) + κ(θ−vt)uv(t,ˆ xt,vt) +1 2vn(t)1−ρ2uxx(t,ˆ xtvt) + 1 2vtσ2uvv(t,ˆ xt,vt)dt =Ztn+1 tn Eρκ σ−1 2(vn(t)−vt)ux(t,ˆ xt,vt) + 1 2(vn(t)−vt)1−ρ2uxx(t,ˆ xt,vt)dt. Since vn(t)−vt=−Zt η(t)κ(θ−vs)ds −σZt η(t)√vsdWs, we have en=e(1) n+e(2) nwith e(1) n:=ρκ σ−1 2Ztn+1 tn E−Zt η(t)κ(θ−vs)ds −σZt η(t)√vsdWsux(t,ˆ xt,vt)dt, e(2) n:=1 21−ρ2Ztn+1 tn E−Zt η(t)κ(θ−vs)ds −σZt η(t)√vsdWsuxx(t,ˆ xt,vt)dt. By Theorem 3and Corollary 1, we have that ux and uxx are bounded. So, Lemma 1 implies that Ztn+1 tn EZt η(t)κ(θ−vs)ds ux(t,ˆ xt,vt)dt =O((∆t)2) and Ztn+1 tn EZt η(t)κ(θ−vs)ds uxx(t,ˆ xt,vt)dt =O((∆t)2).
Risks 2021,9, 23 31 of 38 5.1. Euler Scheme: Preliminaries By the Lemmata 1,4, and 7, we have that sup t∈[0,T] E|vt|p+sup t∈[0,T] E|ˆ xt|p<∞ and E|vt−vs|p+E|ˆ xt−ˆ xs|p≤c·|t−s|p/2,s,t∈[0, T], for all p≥ 1. Using the Burkholder-Davis-Gundy, Hölder, and Minkowski inequalities, we also have sup t∈[0,T] E|Xt|p<∞ and E|Xt−Xs|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) n:=Ztn+1 tn Eρκ σ−1 2−Zt η(t)κ(θ−vs)ds −σZt η(t)√vsdWsux(t,ˆ xt,vt)dt, (19) e(2) n:=Ztn+1 tn E1−ρ2 2−Zt η(t)κ(θ−vs)ds −σZt η(t)√vsdWsuxx(t,ˆ xt,vt)dt, (20) in Section 4.1. If higher derivatives of uare available, then we can analyze EZt η(t)√vsdWsux(t,ˆ xt,vt) and EZt η(t)√vsdWsuxx(t,ˆ xt,vt) via another application of It ¯ o’s lemma. So, let k:[ 0, T]×R×[ 0, ∞)→R be a C1,2 -function that fulfills the backward PDE (7) . In particular, the partial derivatives of u up to order two are such functions. It¯ o’s formula and the Kolmogorov backward PDE (7) now give k(t,ˆ xt,vt) = k(η(t),ˆ xη(t),vη(t)) +Zt η(t)"ρκ σ−1 2kx(s,ˆ xs,vs) + 1−ρ2 2kxx(s,ˆ xs,vs)#(vη(s)−vs)ds +Zt η(t)kx(s,ˆ xs,vs)q1−ρ2vη(s)dBs+Zt η(t)kv(s,ˆ xs,vs)σ√vsdWs. If kx and kv have polynomial growth, then an application of the It ¯ o isometry and the martingale property of the It¯ o integral yield
Risks 2021,9, 23 32 of 38 EZt η(t)√vsdWsk(t,ˆ xt,vt)=EZt η(t)√vsdWsk(η(t),ˆ xη(t),vη(t)) +EZt η(t)√vsdWsZt η(t)K(s,ˆ xs,vs)(vη(s)−vs)ds +EZt η(t)√vsdWsZt η(t)kx(s,ˆ xs,vs)q1−ρ2vη(s)dBs +EZt η(t)√vsdWsZt η(t)kv(s,ˆ xs,vs)σ√vsdWs =EZt η(t)√vsdWsZt η(t)Zs η(s)K(s,ˆ xs,vs)κ(vu−θ)duds −σEZt η(t)√vsdWsZt η(t)K(s,ˆ xs,vs)Zs η(s)√vudWuds +σEZt η(t)vskv(s,ˆ xs,vs)ds, where K(s,ˆ xs,vs) = ρκ σ−1 2kx(s,ˆ xs,vs) + 1−ρ2 2kxx(s,ˆ xs,vs). If kx and kxx have polynomial growth, then an application of Hölder’s inequality and the It¯ o isometry yield EZt η(t)√vsdWsZt η(t)Zs η(s)K(s,ˆ xs,vs)κ(vu−θ)duds=O((∆t)5/2), and so it follows EZt η(t)√vsdWsk(t,ˆ xt,vt)=σEZt η(t)vskv(s,ˆ xs,vs)ds −σEZt η(t)K(s,ˆ xs,vs)Zt η(t)√vudWuZs η(s)√vudWuds +O((∆t)5/2). Since we have EZt η(t)K(s,ˆ xs,vs)Zt η(t)√vudWuZs η(s)√vudWuds =EZt η(t)K(s,ˆ xs,vs)EZt η(t)√vudWuZs η(s)√vudWuFsds =E"Zt η(t)K(s,ˆ xs,vs)Zs η(s)√vudWu2 ds#, again, by the properties of the It ¯ o integral, we finally obtain by Hölder’s inequality and the Burkholder-Davis-Gundy inequality that EZt η(t)K(s,ˆ xs,vs)Zt η(t)√vudWuZs η(s)√vudWuds=O((∆t)2). Thus, we can conclude that EZt η(t)√vsdWsk(t,ˆ xt,vt)=σEZt η(t)vskv(s,ˆ xs,vs)ds+O((∆t)2), (21) for k=uxand k=uxx, if the derivatives up to order four of uhave polynomial growth.
Risks 2021,9, 23 33 of 38 5.2. Euler Scheme: Conclusion Setting now k=uxin Equation (21), we have from (19) that e(1) n:=−ρκ σ−1 2Ztn+1 tnZt η(t) E[ux(t,ˆ xt,vt)κ(θ−vs)]dsdt −σ2ρκ σ−1 2Ztn+1 tnZt η(t) E[vsuxv(s,ˆ xs,vs)]dsdt +O((∆t)3). Replacing now the function kby uxx in Equation (21), we arrive from (20) at e(2) n:=−1 21−ρ2Ztn+1 tnZt η(t) E[uxx(t,ˆ xt,vt)κ(θ−vs)]dsdt −σ2 21−ρ2Ztn+1 tnZt η(t) E[vsuxxv(s,ˆ xs,vs)]dsdt +O((∆t)3). Summarizing, we have shown that en=1 2−ρκ σZtn+1 tn EZt η(t)κ(θ−vs)ux(t,ˆ xt,vt) + σ2vsuxv(s,ˆ xs,vs)dsdt −(1−ρ2) 2Ztn+1 tn EZt η(t)κ(θ−vs)uxx(t,ˆ xt,vt) + σ2vsuxxv(s,ˆ xs,vs)dsdt +O((∆t)3) and E[f(xN,vN)]−E[f(XT,VT)]= N−1 ∑ n=0Ztn+1 tnZt tn EH(s,t,ˆ xs,ˆ xt,vs,vt)dsdt +O((∆t)2) where H(s,t,ˆ xs,ˆ xt,vs,vt) = 1 2−ρκ σκ(θ−vs)ux(t,ˆ xt,vt) + σ2vsuxv(s,ˆ xs,vs) −(1−ρ2) 2κ(θ−vs)uxx(t,ˆ xt,vt) + σ2vsuxxv(s,ˆ xs,vs). 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,6yields for s,t∈[tn,tn+1]that E[H(s,t,ˆ xs,ˆ xt,vs,vt)]=E[H(tn,tn,Xtn,Xtn,Vtn,Vtn)]+O((∆t)1/4). Note here that u∈C4 pol,T(R×R+;R) implies that utx and utxx are well-defined, have polynomial growth, and are continuous. Thus, for an equidistant discretization tk=kT/N,k=0, . . . , N, we have E[f(xN,vN)]−E[f(XT,VT)]=∆t 2 N−1 ∑ n=0 E[H(tn,tn,Xtn,Xtn,Vtn,Vtn)]∆t+O((∆t)5/4). Since N−1 ∑ n=0 E[H(tn,tn,Xtn,Xtn,Vtn,Vtn)]∆t→ZT 0 E[H(t,t,Xt,Xt,Vt,Vt)]dt for ∆t→0, this concludes the proof of Theorem 2(i).
Risks 2021,9, 23 34 of 38 5.3. Semi-Trapezoidal Scheme: Preliminaries By the Lemmata 1,8, and 11, we have that sup t∈[0,T] E|vt|p+sup t∈[0,T] E|ˆ xt|p<∞ and E|vt−vs|p+E|ˆ xt−ˆ xs|p≤c·|t−s|p/2,s,t∈[0, T], for all p≥ 1. Using the Burkholder-Davis-Gundy, Hölder, and Minkowski inequalities, we also have sup t∈[0,T] E|Xt|p<∞ and E|Xt−Xs|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. We will now take also a closer look at the error of the semi-trapezoidal discretization for u∈C4 pol,T(R×R+;R). Recall that e(1) n:=ρκ σ−1 2Ztn+1 tn E1 2−Zt η(t)κ(vt−vs)ds −σZt η(t)√vsdWsux(t,ˆ xt,vt), e(2) n:=Ztn+1 tn E1 2(1−ρ2)−Zt η(t)κ(θ−vs)ds −σZt η(t)√vsdWsuxx(t,ˆ xt,vt) +ρκ σ−1 22σ2 8(t−η(t))2vsuxx(t,ˆ xt,vt)#dt, e(3) n:=ρκ σ−1 2σ2 2Ztn+1 tn E[(t−η(t))vtuxv(t,ˆ xt,vt)]dt. We can again use the It ¯ o formula and the Kolmogorov backward PDE (7) evaluated at (s,ˆ xs,vs)and obtain for a C1,2-function k, which fulfills the PDE (7), that k(t,ˆ xt,vt)−k(tn,ˆ xtn,vtn) =Zt tn kt(s,ˆ xs,vs)ds +Zt tn kv(s,ˆ xs,vs)dvs+Zt tn kx(s,ˆ xs,vs)dˆ xs +1 2Zt tn kxx(s,ˆ xs,vs)dhˆ xis+Zt tn kxv(s,ˆ xs,vs)dhˆ x,vis+1 2Zt tn kvv(s,ˆ xs,vs)dhvis =Zt tnas−ρκ σ−1 2vskx(s,ˆ xs,vs)ds +1 2Zt tnb2 s+c2 s−(1−ρ2)vskxx(s,ˆ xs,vs)ds +Zt tn csσ√vskxv(s,ˆ xs,vs)ds +Zt tn bskx(s,ˆ xs,vs)dBs+Zt tn cskx(s,ˆ xs,vs)dWs +Zt tn σ√vskv(s,ˆ xs,vs)dWs (22) with at:=ρκ σ−1 2vη(t)+1 2(t−η(t))κ(θ−vt) + 1 2(vt−vη(t)), bt:=q1−ρ2pvη(t), ct:=ρκ σ−1 21 2(t−η(t))σ√vt.
Risks 2021,9, 23 35 of 38 Analogous calculations as for the Euler scheme yield that EZt η(t)√vτdWτk(t,ˆ xt,vt)=σEZt η(t)vτkv(τ,ˆ xτ,vτ)dτ+O((∆t)2)(23) for k=uxand k=uxx under the assumption u∈C4 pol,T(R×R+;R). 5.4. Semi-Trapezoidal Rule: Calculations for e(1) n, e(2) n, and e(3) n Rewriting the terms of e(1) nusing (23) for the last term gives e(1) n=−ρκ σ−1 21 2Ztn+1 tn EZt η(t)κ(vt−vs)ds ux(t,ˆ xt,vt)dt −ρκ σ−1 2σ 2Ztn+1 tn EZt η(t)√vsdWsux(t,ˆ xt,vt)dt =−ρκ σ−1 21 2Ztn+1 tn EZt η(t)κZt sκ(θ−vu)du +σZt s√vudWuds ux(t,ˆ xt,vt)dt −ρκ σ−1 2σ2 2Ztn+1 tnZt η(t) E[vsuxv(s,ˆ xs,vs)]dsdt +O((∆t)3). Applying, again, (23) with s instead of η(t) 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 Ztn+1 tn EZt η(t)κZt sκ(θ−vu)du +σZt s√vudWuds ux(t,ˆ xt,vt)dt =O((∆t)3) and so e(1) n=−ρκ σ−1 2σ2 2Ztn+1 tnZt η(t) E[vsuxv(s,ˆ xs,vs)]dsdt +O((∆t)3). Adding e(3) nyields e(1) n+e(3) n=O((∆t)3)−ρκ σ−1 2σ2 2Ztn+1 tnZt tn E[uxv(t,ˆ xt,vt)vt−uxv(s,ˆ xs,vs)vs]dsdt. However, It¯ o’s formula gives for sufficiently smooth k:[0, T]×R×[0, ∞)→Rthat k(t,ˆ xt,vt)−k(s,ˆ xs,vs) =Zt skt(r,ˆ xr,vr)dr +Zt skv(r,ˆ xr,vr)dvr+Zt skx(r,ˆ xr,vr)dˆ xr +1 2Zt skxx(r,ˆ xr,vr)dhˆ xir+Zt skxv(r,ˆ xr,vr)dhˆ x,vir+1 2Zt tn kvv(r,ˆ xr,vr)dhvir (24) with dvt=κ(θ−vt)dt +σ√vtdWt,dˆ xt=atdt +btdBt+ctdWt,
Risks 2021,9, 23 36 of 38 where at:=ρκ σ−1 2vη(t)+1 2(t−η(t))κ(θ−vt) + 1 2(vt−vη(t)), bt:=q1−ρ2pvη(t), ct:=ρκ σ−1 21 2(t−η(t))σ√vt and dhˆ xit= (b2 t+c2 t)dt,dhˆ x,vit=σct√vtdt,dhvit=σ2vtdt. Since u∈C4 pol,T(R×R+;R) , we can apply this to k(t , x , v) = uxv(t , x , v)v and taking expectations gives then E[uxv(t,ˆ xt,vt)vt−uxv(s,ˆ xs,vs)vs]=O(|t−s|). So, we end up with e(1) n+e(3) n=O((∆t)3). Looking at e(2) n, the last term is already of third order: e(2) n=Ztn+1 tn E1 2(1−ρ2)−Zt η(t)κ(θ−vs)ds −σZt η(t)√vsdWsuxx(t,ˆ xt,vt) +ρκ σ−1 22σ2 8(t−η(t))2vsuxx(t,ˆ xt,vt)#dt =Ztn+1 tn E1 2(1−ρ2)−Zt η(t)κ(θ−vs)ds −σZt η(t)√vsdWsuxx(t,ˆ xt,vt)dt +O((∆t)3). Since EZt η(t)√vsdWsuxx(t,ˆ xt,vt)=σEZt η(t)vsuxxv(s,ˆ xs,vs)ds+O((∆t)2), by (23), it follows e(2) n=−1 2(1−ρ2)Ztn+1 tn EZt η(t)κ(θ−vs)ds uxx(t,ˆ xt,vt)dt −1 2(1−ρ2)σ2Ztn+1 tnZt η(t) E[vsuxxv(s,ˆ xs,vs)]dsdt +O((∆t)3). 5.5. Semi-Trapezoidal Scheme: Conclusion Summarizing, we have shown that en=−(1−ρ2) 2Ztn+1 tn EZt η(t)κ(θ−vs)uxx(t,ˆ xt,vt) + σ2vsuxxv(s,ˆ xs,vs)dsdt +O((∆t)3) and E[f(xN,vN)]−E[f(XT,VT)]= N−1 ∑ n=0Ztn+1 tnZt tn EH(s,t,ˆ xs,ˆ xt,vs,vt)dsdt +O((∆t)2)
Risks 2021,9, 23 37 of 38 where H(s,t,ˆ xs,ˆ xt,vs,vt) = −(1−ρ2) 2κ(θ−vs)uxx(t,ˆ xt,vt) + σ2vsuxxv(s,ˆ xs,vs). 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∈[tn,tn+1]that E[H(s,t,ˆ xs,ˆ xt,vs,vt)]=E[H(tn,tn,Xtn,Xtn,Vtn,Vtn)]+O((∆t)1/4). In particular, for an equidistant discretization tk=kT/N,k=0, . . . , N, we have E[f(xN,vN)]−E[f(XT,VT)]=∆t 2 N−1 ∑ n=0 EH(tn,tn,Xtn,Xtn,Vtn,Vtn)∆t+O((∆t)5/4) and the convergence N−1 ∑ n=0 EH(tn,tn,Xtn,Xtn,Vtn,Vtn)∆t→ZT 0 EH(t,t,Xt,Xt,Vt,Vt)dt for ∆t→0 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. [CrossRef] Altmayer, Martin. 2015. Quadrature of Discontinuous SDE Functionals Using Malliavin Integration by Parts. Ph.D. dissertation, University of Mannheim, Germany. 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. 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. [CrossRef] Andersen, Leif. 2008. Simple and efficient simulation of the Heston stochastic volatility model. Journal of Computational Finance 11: 29–50. 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. [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. 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. Broadie, Mark, and Özgür Kaya. 2006. Exact simulation of stochastic volatility and other affine jump diffusion processes. Operations Research 54: 217–31. [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. 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. 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. [CrossRef]
Risks 2021,9, 23 38 of 38 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. [CrossRef] Glasserman, Paul, and Kyoung-Kuk Kim. 2011. Gamma expansion of the Heston stochastic volatility model. Finance and Stochastics 15: 267–96. [CrossRef] 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. [CrossRef] Hurd, Thomas, and Alexey Kuznetsov. 2008. Explicit formulas for Laplace transforms of stochastic integrals. Markov Processes and Related Fields 14: 277–90. Kahl, Christian, and Peter Jäckel. 2006. Fast strong approximation Monte-Carlo schemes for stochastic volatility models. Quantitative Finance 6: 513–36. Karatzas, Ioannis, and Steven Shreve. 1991. Brownian Motion and Stochastic Calculus, 2nd ed. New York: Springer. 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. Nualart, David. 1995. The Malliavin Calculus and Related Topics. New York: Springer. Lord, Roger, Remmert Koekkoek, and Dick van Dijk. 2009. A comparison of biased simulation schemes for stochastic volatility models. Quantitative Finance 10: 177–94. Smith, Robert. 2007. An almost exact simulation method for the Heston model. Journal of Computational Finance 11: 115–25. 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. Zheng, Chao. 2017. Weak convergence rate of a time-discrete scheme for the Heston stochastic volatility. SIAM Journal on Numerical Analysis 55: 1243–63. [CrossRef]