1 Introduction
In this paper, we study the total variation distance between a one-dimensional diffusion process and a third weak order scheme based on a cubature method and introduced by Alfonsi [1]. In his work, Alfonsi proved that it converges with weak order 3 for smooth test functions with polynomial growth. We will show that the convergence also takes place with order 3 if we consider measurable and bounded test functions. In this case, we say that the total variation distance between the diffusion process and the scheme converges towards zero with order 3. In order to do it, we will use a result from [3] based on en abstract Malliavin calculus introduced by Bally and Clément [2]. A main interest of this approach is that the random variables used to build the scheme are not necessarily Gaussian but belong to a class of random variables with no specific law. Consequently, our result can be seen as an invariance principle.
Let us be more specific. We consider the ℝ-valued one-dimensional Markov diffusion process
(1.1)dXt=V0(Xt)dt+V1(Xt)∘dWt,
where Vi∈𝒞b∞(ℝ,ℝ) for i=0,1, (Wt)t⩾0 is a one-dimensional standard Brownian motion and ∘dWt is the Stratonovich integral with respect to Wt. In this paper, we will study an approximation scheme for (1.1) which is defined on an homogeneous time grid. It is relevant to notice that the results we will obtain remain true for non-homogeneous time grids, but we do not treat that case for sake of clarity. We fix T>0 and we denote n∈ℕ*, the number of time step between 0 and T. Then, for k∈ℕ we define tkn=kT/n and we introduce the homogeneous time grid πT,n={tkn=kT/n:k∈ℕ} and its bounded version πT,nT~={t∈πT,n:t⩽T~} for T~⩾0. Finally, for S∈[0,T~) we will denote πT,nS,T~={t∈πT,nT~:t>S}. Now, for tkn=kT/n, we introduce the abstract ℝ-valued Markov chain
(1.2)Xtk+1nn=ψk(Xtknn,Zk+1n,δk+1n),k∈ℕ,
where ψk:ℝ×ℝ×ℝ+→ℝ is a smooth function such
that ψk(x,0,0)=x, Zk+1∈ℝN, k∈ℕ, is a sequence of independent and centered
random variables and supk∈ℕ∗δkn⩽C/n.
Before estimating the distance between X and Xn, we introduce some notations. For f∈𝒞∞(ℝd) and for a multi-index α=(α1,…,αd)∈ℕd we denote |α|=α1+…+αd and ∂αf=∂xαf=∂x1α1…∂xdαdf(x). We include
the multi-index α=(0,…,0) and in this case ∂αf=f. We will use the norms
∥f∥q,∞=supx∈ℝd∑0⩽|α|⩽q|∂αf(x)|,q∈ℕ.
In particular, ∥f∥0,∞=∥f∥∞ is the usual supremum norm and we denote 𝒞bq(ℝd)={f∈𝒞q(ℝd):∥f∥q,∞<∞}.
A first standard result is the following: Let us assume that there exist h>0 and q∈ℕ such
that for every test function f∈𝒞bq(ℝ), k∈ℕ and x∈ℝ,
(1.3)|𝔼[f(Xtk+1n)-f(Xtk+1nn)|Xtkn=Xtknn=x]|⩽C∥f∥q,∞nh+1.
Then we have
(1.4)supt∈πT,nT|𝔼[f(Xt)-f(Xtn)]|⩽C∥f∥q,∞nh.
It means that (Xtkn)k∈ℕ is an approximation scheme of weak order h for the Markov process (Xt)t⩾0 for the test functions f∈𝒞bq(ℝ;ℝ). The value h thus measures the efficiency of the scheme whereas q stands for the required regularity on the test functions in order to obtain convergence with order h. This subject has already been studied widely in the literature and we point out some famous examples. However, the reader may notice that in all those works, the required order of regularity q is greater than one. Concerning the Euler scheme for diffusion processes, the result (1.4), with h=1, has initially been proved in the seminal papers of Milstein [19] and of
Talay and Tubaro [22] (see also [12]). Since then, various situations have been studied: Diffusion
processes with jumps (see [21, 10]) or diffusion processes with boundary
conditions (see [7, 6, 8]). An overview of the subject is proposed in [11].
More recently, discretization schemes of higher orders (e.g., h=2), based on
cubature methods, have been introduced and studied by Kusuoka [16], Lyons
[18], Ninomiya and Victoir [20] or Alfonsi [1]. The reader may also refer to the work Kohatsu-Higa and Tankov [13] for a higher weak order for jump processes. Finally, in [1], a third weak order scheme (with h=3) has been introduced following similar cubature ideas. This is the one we will study in this paper.
As we already made precise, all those schemes converge for some q⩾1 in (1.4). Another point of interest relies thus on the study of the set of test functions which enable the converge with weak order h. The purpose is to extend this set beyond 𝒞bq(ℝ;ℝ) and to obtain (1.4) with ∥f∥q,∞ replaced by ∥f∥∞ when f is a measurable and bounded function. In this case, we say that the scheme converges for the total variation distance. A first result of this type has been obtained
by Bally and Talay [4, 5]. They treat the case of the Euler scheme using the Malliavin calculus (see also
Guyon [9] when f is a tempered distribution). Afterwards Konakov, Menozzi and Molchanov [14, 15]
established some local limit theorems using a parametrix method. Recently Kusuoka [17], also using Malliavin Calculus,
obtained estimates of the error in total variation distance for the Ninomiya Victoir scheme (which corresponds to the case h=2) under a Hörmander-type condition.
Under an ellipticity condition, we will obtain a similar result for the case h=3, using a scheme introduced in [1]. This scheme is well defined if the Lie bracket between V12 and V0 is equal to 2V~2, with V~ a first order differential operator. Since we consider one-dimensional processes with form (1.1), we will be able to give an explicit necessary and sufficient condition in order to obtain this property on the Lie bracket.
Moreover, we will not work in a Gaussian framework and then we will have to use a variant of the Malliavin calculus introduced by Bally and Clement [2] for which we can apply the results from [3]. A main interest of this approach is that the random variables involved in the scheme do not have a specific law but simply belong to a class of random variables which are Lebesgue lower bounded and satisfy some moment conditions. In this way, our final result can be seen as an invariance principle. The ambit of this scheme thus goes well beyond the Gaussian case.
We will begin presenting the framework of this paper in Section 2. In Section 3, we will give some third weak order convergence results for smooth test functions and for bounded measurable test functions. The latter is presented in Theorem 3.2 and constitutes the main result of this paper. It gives the convergence for the total variation distance with order 3 of the scheme from [1], toward the Markov process (1.1). We will also obtain an estimate of the density function of the diffusion and its derivatives. We will follow with a short numerical illustration in order to check the order of convergence for a suited example. This paper will end with the proof of our main theorems in Section 6.
2 The third weak order scheme
We consider the one-dimensional ℝ-valued diffusion process
(2.1)dXt=V0(Xt)dt+V1(Xt)∘dWt
with V0,V1∈𝒞b∞(ℝ;ℝ), (Wt)t⩾0 a standard Brownian motion. Moreover, ∘dWt denotes the
Stratonovich integral with respect to W. The infinitesimal operator of
this Markov process is
A=V0+12V12
with the notation Vf(x)=V(x)∂f(x). Let us define exp(V)(x):=ΦV(x,1), where ΦV solves the deterministic equation
(2.2)ΦV(x,t)=x+∫0tV(ΦV(x,s))𝑑s.
By a change of variables one obtains ΦεV(x,t)=ΦV(x,εt), so we have
exp(εV)(x):=ΦεV(x,1)=ΦV(x,ε).
We also notice that the semigroup of the above Markov process, which is given by PtVf(x)=f(ΦV(x,t)), has the infinitesimal operator AVf(x)=Vf(x). In particular, the relation PtVAV=AVPtV
reads
Vf(ΦV(x,t))=AVPtVf=PtVAVf=V(x)∂x(f∘ΦV)(x,t).
Using m times Dynkin’s formula PtVf(x)=f(x)+∫0tPsVAVf(x)𝑑s, we obtain
(2.3)f(ΦV(x,t)))=f(x)+∑r=1mtrr!Vrf(x)+1m!∫0t(t-s)mVm+1PsVf(x)ds.
We present now the third weak order scheme introduced in [1]. In order to do it, we introduce the following commutation property:
(2.4)V12V0-V0V12=2V~2,
where V~ is a first order operator. We consider some sequences ϵk,ρk, k∈ℕ, of independent uniform random variables with values in {-1,1} and {1,2,3}, and we define ψ:{-1,1}×{1,2,3}×ℝ3→ℝ using the following splitting procedure:
(2.5)ψ(ϵk,ρk,x,wk+11,wk+10)={exp(ϵkwk+10V~)∘exp(wk+10V0)∘exp(wk+11V1)(x),if ρk=1,exp(wk+10V0)∘exp(ϵkwk+10V~)∘exp(wk+11V1)(x),if ρk=2,exp(wk+10V0)∘exp(wk+11V1)∘exp(ϵkwk+10V~)(x),if ρk=3,
with wk0=T/n, wk1=TZk/n. We notice that ψ(ϵk,ρk,x,0,0)=x, which is relevant with the definition of a scheme.
Moreover, Zk, k∈ℕ∗, are independent random variables which
are lower bounded by the Lebesgue measure: There exist z∗,k∈ℝ and ε∗,r∗>0 such that for every Borel set A⊂ℝ and every k∈ℕ∗,
(2.6)Lz∗(ε∗,r∗) ℙ(Zk∈A)⩾ε∗λ(A∩Br∗(z∗,k)).
Moreover, we assume that the sequence Zk satisfies the following moment conditions:
(2.7)𝔼[Zk]=𝔼[Zk3]=𝔼[Zk5]=𝔼[Zk7]=0, 𝔼[Zk2]=1,𝔼[Zk4]=3,𝔼[Zk6]=15,𝔼[|Zk|p]<∞,p⩾1.
One step of our scheme (between times tkn and tk+1n) is given by
(2.8)Xtk+1nn=ψ(ϵk,ρk,Xtknn,wk+11,wk+10).
Using the notation from (1.2), we also have
(2.9)Xtk+1nn=ψk(Xtknn,wk+11,wk+10)
with ψk(x,z,t)=ψ(ϵk,ρk,x,z,t). In the sequel, we will study the third order convergence of this scheme towards the Markov process given in (1.1) for smooth test functions and for bounded measurable test functions.
3 Convergence results
We begin introducing some notations. Let r∈ℕ∗. For a sequence of functions ψk∈𝒞r(ℝ×ℝ×ℝ+;ℝ), k∈ℕ, we denote
∥ψ∥1,r,∞=1 and supk∈ℕ∑|α|=0r∑|β|+|γ|=1r-|α|∥∂xα∂zβ∂tγψk∥∞,
and for r∈ℕ∗,
(3.1)𝔎r(ψ)=(1+∥ψ∥1,r,∞)exp(∥ψ∥1,3,∞2).
3.1 Smooth test functions
In this section, we study the convergence of the scheme given in (2.9) for smooth test functions. We state a first result, which is the starting point in order to prove the convergence in total variation distance.
Theorem 3.1Suppose that V0,V1,V~∈Cb∞(R;R). We also assume that (2.4) and (2.7) hold. Then there exists some universal constants l∈N∗ and C⩾1
such that for every f∈Cb8(R), we have
(3.2)supt∈πT,nT|𝔼[f(Xt))-𝔼[f(Xtn)]|⩽CC8(V)l∥f∥8,∞n3
with Cq(V):=supi=0,1∥Vi∥q,∞+∥V~∥q,∞.
Remark 3.1This result has already been obtained in [1] in the case of test functions with polynomial growth. The proof is similar and since we intend to obtain this result with the supremum norm of f we do not treat that case.
We give a proof of this result in Section 6. Once we have used the Lindeberg decomposition, it relies on short time estimates using the Dynkin’s formula. Now, we are going to take a step further and consider simply bounded and measurable test functions. Notice that, it means the convergence for total variation distance.
3.2 Bounded measurable test functions
We see that estimate (3.2) involves the derivatives of order 8 of the test function. We will see that it is possible to obtain similar estimates with ∥f∥8,∞ replaced by ∥f∥∞. This is a consequence of a result from [3] in which the authors provide some sufficient conditions for the scheme in order to obtain the convergence for the total variation distance. The scheme (2.8) satisfy those conditions and, under an ellipticity assumption on the diffusion coefficient V1, we are going to obtain an estimate of its total variation distance with the diffusion process (2.1).
Before doing it, we introduce a necessary and sufficient explicit condition in order to obtain (2.4) as soon as for all x∈ℝ,V1(x)≠0. Notice that, since we assume that V1 is continuous, it has a constant sign. Moreover, this hypothesis will not be restrictive in this application. Indeed, the ellipticity condition required to use the result from [3] implies that infxV1(x)2⩾λ∗>0 for a constant λ∗. We will suppose without loss of generality that V1 is positive. The necessary and sufficient condition for (2.4) is the following: We assume that the function
(3.3)g:ℝ→ℝ,x↦V0(x)/V1(x)
is increasing. Notice that if V1 is negative, g has to be decreasing.
Moreover, we propose an alternative scheme in order to approximate the density function of X and its derivatives. We consider a standard normal random variable
G which is independent from Zk,k∈ℕ, and for θ>0, we introduce (Xtn,θ)t∈πT,n as follows:
Xtn,θ(x)=1nθG+Xtn(x),
where Xn(x) is the process which starts from x that is X0n=x. We denote by ptθ,n(x,y) the density of the law of Xtθ,n(x) and for t∈πT,n, we define
(3.4)Qtn,θf(x):=𝔼[f(1nθG+Xtn(x))].
Now, we can state our main result.
Theorem 3.2Suppose that V0,V1,V~∈Cb∞(R;R). We fix T>0 and we also assume that (3.3), (2.6) and (2.7) hold and that
(3.5)V1(x)2⩾λ∗>0 for all x∈ℝ.
Let S∈(0,T/2). Then there exists n0∈N∗ such that for every n⩾n0, we have the following properties.
There exist l∈ℕ∗ and C⩾1 which depends on m∗,
r∗ and the moments of Z such that, for every bounded and measurable function
f:ℝ→ℝ,
(3.6)supt∈πT,n2S,T|𝔼[f(Xt)]-𝔼[f(Xtn)]|⩽CC8(V)l𝔎11(ψ)l(λ∗S)η(8)∥f∥∞n3
with 𝔎r(ψ) and Cq(V) given in (3.1) and (3.2) and η(r)=r(r+1).
Moreover, for every t>0,Pt(x,dy)=pt(x,y)dy with (x,y)↦pt(x,y) belonging to 𝒞∞(ℝd×ℝd).
Let θ⩾h+1. We recall the Qn,θ is defined in (3.4) and verifies Qtn,θ(x,dy)=ptn,θ(x,y)dy. Then, there exists l∈ℕ∗ such that for every R>0,ε∈(0,1),
x0,y0∈ℝd, and every multi-index α,β with |α|+|β|=u, we also have
supt∈πT,n2S,Tsup(x,y)∈B¯R(x0,y0)|∂xα∂yβpt(x,y)-∂xα∂yβptn,θ(x,y)|⩽CC8(V)l𝔎11(ψ)l(λ∗S)η(pu,ε∨8)1n3(1-ε)
with a constant C which depends on R,
x0,
y0,
T,
|α|+|β| and pu,ε=(u+2d+1+2⌈(1-ε)(u+d)/(2ε)⌉).
Remark 3.2It is relevant to notice that we have the same result if we assume that the function defined in (3.3) is decreasing (resp. increasing) for V1 positive (resp. V1 negative). In this case V0V12-V12V0=2V~2 and we have to define the scheme differently. In the construction (2.5), we invert the terms containing V1 with the ones containing V0.
Remark 3.3Property (2.6) is crucial here, since we will use a result from [3] which employs abstract integration by parts formulae based on the noise Zk. However it is not restrictive for concrete applications.
The result (3.6) signifies the convergence in total variation with order 3. The proof of this theorem is given in Section 6. Since we have already obtained some short time estimates of the form (1.3) in the proof of Theorem 3.1 and (3.3) holds, the key point of this proof does not rely on the weak order of the scheme. This is the fact that, the splitting procedure (2.5) in order to build the scheme, always includes a diffusion part through exp(Zk/n/TV1), with Zk satisfying (2.6) and the ellipticity condition (3.5) for V1. The proof is then a consequence from [3, Theorem 3.3] which employs an abstract Malliavin calculus based on such noise Zk and initially presented by Bally and Clément [2]. A similar approach can be used in order to prove the convergence for the total variation distance for even higher order scheme built as in (2.5). The main difficulty will then rely on the proof of the short time estimate (1.3).
A main interest of this result is that it can be seen as an invariance principle as well. Indeed, it does not require that Zk follows a particular law but only properties (2.6) and (2.7). In particular, we do not restrict ourselves to the Gaussian framework which is necessary to use the Malliavin Calculus in order to prove the convergence for the total variation distance as in [4], [5], or [17]. In this way, condition (2.6) might be a hint to find a necessary condition on the random variables (Zk)k∈ℕ* in order to obtain the total variation convergence with order h=3.
Moreover, using Remark 3.2, we can define a third order scheme as soon as the function defined in (3.3) is monotonic. If it is increasing (recall that V1(x)⩾0), the Lie bracket between V12 and V0 is given by [V12,V0]f=V12V0f-V0V12f=2V~2f with
V~(x)=|V1(x)(V1(x)∂xV0(x)-∂xV1(x)V0(x))|,x∈ℝ.
If it is decreasing, we have [V0,V12]f=2V~2f as well. This explicit representation for V~ is crucial for concrete applications since the scheme is defined using the solution of (2.2) with V=V~. Moreover, looking at (3.2) and (3.6), we have to control its derivatives.
4 Numerical illustration
In this section, we study the numerical approximation of a one-dimensional SDE with schemes defined on homogeneous time grids with form πT,n={kT/n:k∈ℕ}. We will fix T and we will analyze the behavior of the total variation distance between the diffusion process (Xt)t⩾0 and miscellaneous discretization schemes (Xtn)t∈πT,n with respect to the number of time step n. More particularly, we will study the weak error |𝔼[f(Xt)]-𝔼[f(Xtn)]| for bounded measurable functions f and various n.
In concrete applications, once we have selected a scheme Xn, 𝔼[f(Xtn)] will be used to estimate 𝔼[f(Xt)]. The next step is thus to approximate 𝔼[f(Xtn)]. A standard way to do it, is to use a Monte Carlo method. Given an independent sampling of size M, and using the Central Limit Theorem, we can easily show that those algorithms converge toward the real expectancy with rate M. Moreover, discretization schemes provide an estimation of 𝔼[f(Xt)] with any desired precision since we can choose any value for n. However, the cost of calculation will also increase with n since we have n iterations of the scheme function (1.2). At this point, it is important to notice that there is a trade off to make between the precision we want to obtain and the time of calculation we can afford. Indeed, if our scheme converges with order h, we have to choose M=𝒪(n2h) and then choose n large enough in order to obtain the desired precision. We will see that even if the time of calculation of one step of the scheme we study in this paper is much longer than the time of a lower order scheme (e.g., the Euler scheme), the third weak order scheme is better in time of calculation and precision as soon as the precision is high enough. In order to illustrate the reason why we point out such properties, we now present our example.
We consider the Markov diffusion process (Xt)t⩾0 given by the following SDE:
(4.1)dXt=adt+σarctan(Xt)+π∘dWt
with σ>0 and a∈ℝ. Notice that the coefficients of the SDE (4.1) belong to 𝒞b∞(ℝ) and moreover the function V1 given by x↦σ/(arctan(x)+π) satisfies infxV1(x)>2σ/π and the function V0/V1=a/V1 is increasing. Therefore, the scheme (2.5) is well defined and we have the required hypothesis in order to obtain the results from Theorem 3.2. Moreover, we have an explicit representation for the first order operator V~, that is,
V~(x)=σa/(1+x2(arctan(x)+π)3/2).
The next step consists then in solving the ODE (2.2) for V=V0,V1,V~. Looking closer to (2.5), we will use each of these solutions once for each step of the discretization algorithm. In this example, it is easy to find an analytic solution to (2.2) when V=V0. However for V=V1,V~, it is much more cumbersome and we will use some numerical algorithms. A naive algorithm consists in using the Riemann approximation of ∫0tV(ΦV(x,s))𝑑s on a time grid of [0,t] in the following way: For a number N of time steps, we put ΦVN(x,0)=x and for i∈{0,…,N-1}, ΦVN(x,(i+1)t/N)=ΦVN(x,it/N)+TN-1V(ΦVN(x,it/N)). This is the method we will use to approximate ΦV~N(x,t). Finally, in this case, we can use an alternative way to approximate ΦV1N(x,t). Indeed, we can show that g(ΦV1N(x,t))=g(x)+t, where g is the bijective function in 𝒞1(ℝ) defined by
g(x)=xarctan(x)-0.5log(1+x2)+xπσ.
Then, we can find ΦV1(x,t) using a Newton algorithm in order to invert g. Likewise the naive Riemann approximation, this method provides an approximation given a parameter of precision (which is N for Riemann sums). Obviously, the more this parameter is tight, the more the cost of the algorithm is high. Compared to one step the Euler scheme,
Xtk+1nn,Eul=Xtknn,Eul+(a-σ22(1+(Xtknn,Eul)2)(arctan(Xtknn,Eul)+π)3)Tn+σarctan(Xtknn,Eul)+πTNZk+1
with (Zk)k∈ℕ∗ i.i.d ∼𝒩(0,1),
the cost of one step of (2.5) can thus be very important. However, despite that cost, the third order scheme become more effective as soon as we want to compute 𝔼[f(Xt)] with a sufficiently high precision.
Heuristically, let ϵ>0 be the precision of the weak error, that is |𝔼[f(Xt)]-𝔼[f(Xtn)]|⩽ϵ. In order to reach that precision, we will have to run M=ϵ-2 Monte Carlo iterations. Now let n∈ℕ such that n3=ϵ-1 (here, for the sake of clarity, we consider that the right-hand side of estimation (3.6) is of the form 1/n3). Then, if we want to reach this precision, we will have to simulate M=ϵ-2 realizations of the third order scheme with time step t/n, or of the Euler scheme with time step t/n3. Now, we assume that the cost in time of calculation of one step of the third order scheme is given by τNV3 and by τEul for the Euler scheme. Then the total cost to reach the precision ϵ will be τNV3nM=τNV3ϵ-2-1/3 for the third order scheme and τEuln3M=τEulϵ-3 for the Euler scheme. Then, as soon as τNV3/τEul⩽ϵ-2/3, the cost of the third order scheme will be lower than the cost of the Euler scheme. Controversially, if τNV3 and τEul are fixed, we can find a precision ϵ0 such that the cost of the three order scheme is lower than the one of the Euler scheme for all ϵ⩽ϵ0.
In Figure 1, we represent the error |𝔼[f(Xt)]-𝔼[f(Xtn)]|, with respect to the number of time steps n, in Log-Log scale, for the third order scheme we study in this paper and when f is a Heaviside function. We observe that the scheme converges with the expected rate, that is h≈2.91. This numerical experiment thus confirms the total variation convergence result from Theorem 3.2. Notice that we have also implemented the Euler scheme and the Ninomiya Victoir scheme of order 2 (see [20]) in order to compare the cost of the different approaches. With the precision parameters we have selected in the algorithms solving (2.2) in order to obtain Figure 1, we have τNV3≈7.8τNV2≈51.9τEul which is quite reasonable given the gain which is made with respect to the number of time steps. In this case, the third order scheme thus become more effective than the Euler scheme as soon as the precision ϵ of the weak error satisfies |𝔼[f(Xt)]-𝔼[f(Xtn)]|⩽ϵ⩽(51.9)-3/2.
5 Conclusion
In this paper, we provide a practical way to simulate a third weak order scheme for one-dimensional stochastic differential equations. We show that this scheme also converges with third order for the total variation distance. Moreover, using the density of the scheme, we obtain an approximation of the density of the diffusion, with almost third order. Finally, since the random variable which are used to build the scheme only satisfies condition (2.6) and the moment condition (2.6) and then do not have specific law, this result can be seen as an invariance principle.
6 Proof of the main theorems
Proof of Theorem 3.1.
We divide the proof into two steps.
Step 1.
We define (Pt,s)t,s∈πT,n;t⩽s by
Pt,tnf(x)=f(x),Ptkn,tkrnf=𝔼[f(Xtrn)|Xtkn=x]for all k⩽r∈ℕ,
Qt,tnf(x)=f(x),Qtkn,tkrnf=𝔼[f(Xtrnn)|Xtknn=x]for all k⩽r∈ℕ,
and we notice that for t,s,u∈πT,n with t⩽s⩽u, then Pt,unf=Pt,snPs,unf.
It follows that
|𝔼[f(Xtmn)]-𝔼[f(Xtmnn)]|⩽∥P0,tmnnf-Q0,tmnnf∥∞
⩽∑k=0m-1∥Ptk+1n,tmnnPtkn,tk+1nnQtk+1nf-Ptk+1n,tmnnQtkn,tk+1nQtknnf∥∞
=∑k=0m-1∥Ptk+1n,tmn(Ptkn,tk+1n-Qtkn,tk+1n)Qtknf∥∞,
∥Ptmnnf-Qtmnnf∥∞⩽∑k=0m-1∥PtknnPtkn,tk+1nnQtk+1n,tmnnf-PtknnQtkn,tk+1nnQtk+1n,tmnnf∥∞
(6.1)=∑k=0m-1∥Ptknn(Ptkn,tk+1n-Qtkn,tk+1n)Qtk+1n,tmnnf∥∞.
We notice that it easy to prove that, for t,s∈πT,n, t⩽s, ∥Pt,sf∥p,∞⩽C∥f∥p,∞ and ∥Qt,sf∥p,∞⩽C∥f∥p,∞.
Step 2.
It remains to show ∥Ptkn,tk+1nf-Qtkn,tk+1nf∥∞⩽C∥f∥8,∞/n4 and using (6.1) the proof will be completed. In order to simplify the notation, we fix T=1 without loss of generality. For ϵ=-1,1, we denote
𝒯0f(x)=f(exp(1nV0)(x)),𝒯1f(x)=f(exp(ZnV1)(x)),𝒯~ϵf(x)=f(exp(ϵnV~)(x)).
Notice that, using the notation introduced in the beginning of this section with V=n-1/2ZV1, we have 𝒯1f(x)=P1n-1/2ZV1f(x). Using (2.3) with t=1 and V=n-1/2ZV1, we obtain
𝒯1f(x)=f(x)+∑r=1mZrnr/21r!V1rf(x)+Zm+1n(m+1)/2Rm+1,1f(x)
with
(6.2)Rm+1,1f(x)=1m!∫01(1-λ)mV1m+1Pλn-1/2ZV1f(x)𝑑λ
and we recall that Pλn-1/2ZV1f(x)=f(exp(λZV1/n)). We have a similar development if we put V=V0/n or V=ϵV~/n in (2.3).
Our aim is to give a development of order 4 (with respect to n) for 𝔼[f(ψk(x,wk+11,wk+10)] (see (6.3) below). We replace each 𝒯∈{𝒯0,𝒯1,𝒯~ϵ}, with an expansion of order m⩽7 given above with Z=Zk+1 for 𝒯1 and m⩽3 for 𝒯=𝒯0,𝒯~. Then we calculate the products of the miscellaneous expansions, each with a well chosen order such that there is no term with factor n-r, r>4, appearing in those products. Moreover, all the terms containing n-4 go in the
remainder. The last step consists in computing the expectancy. We notice that
𝔼[Ptn-1/2ZV1]=PtV12/(2n)
and 𝔼[Zk+1r]=0 for odd r⩽7. Finally, since 𝔼[Zk+12]=1, 𝔼[Zk+14]=6 and 𝔼[Zk+16]=15, the calculus is completed and we obtain
𝔼[f(ψk(x,wk+11,wk+10)]=16∑ϵ=-1,1𝔼[(𝒯~ϵ𝒯0𝒯1+𝒯0𝒯~ϵ𝒯1+𝒯0𝒯1𝒯~ϵ)f(x)]
=f(x)+1n(V0+12V12)f(x)+12n2(V02+14V14+2V012V12+V~2)f(x)
+16n3(18V16+V03+3V014V14+3V0212V12+2V~212V1
+2V0V~2+12V12V~2+V~2V0)f(x)+1n4Rf(x)
=f(x)+1n(V0+12V12)f(x)+12n2(V0+12V12)2f(x)
(6.3)+16n3(V0+12V12)3f(x)+1n4Rf(x).
The remainder R is a sum of terms of the form
C(𝒯~ϵ,αϵ𝒯0,α0𝒯1,α1+𝒯0,α0𝒯~ϵ,αϵ𝒯1,α1+𝒯0,α0𝒯1𝒯~ϵ,αϵ)f(x)
with α=(α0,α1,α2)∈{0,…,4}3, |α|=α0+α1+α2=4, and, using the notation given in (6.2),
𝒯0,k∈{V0k,Rk,0},𝒯~ϵ,k∈{V~k,R~k,ϵ},𝒯1,k ∈{V12k,R2k,1},k=0,…,3,
𝒯0,4=R4,0,𝒯~ϵ,4=R~4,ϵ,𝒯1,4 =𝐑8,1,
with
𝐑8,1=𝔼[Z8R8,1]=∫01(1-λ)7𝔼[Z8V18PλU1f(x)]𝑑λ.
It is easy to check that for every g∈𝒞k+p(ℝ), we have the following property:
∥𝒯0,kg∥p,∞+∥𝒯1,kg∥p,∞+∥𝒯~ϵ,kg∥p,∞⩽CC2k+p(V)l∥g∥2k+p,∞
for some constants l∈ℕ∗, C⩾1. So
(6.4)∥Rf∥∞⩽CC8(V)l∥f∥8,∞.
We turn now to the diffusion process Xt. We have the development
𝔼[f(Xt(x))]=PtAf(x)=f(x)+tAf(x)+t22A2f(x)+t36A3(x)+t44!Rt′f(x)
with
(6.5)Rt′f(x)=t-1∫0tPλAA4f(x)(1-λt)3𝑑λ.
We take t=n-1 and make the
difference between (6.5) and (6.3). All the terms cancel except
for the remainders, so we obtain
(6.6)𝔼[f(Xtk+1n)]-𝔼[f(Xtk+1nn)∣Xtkn=Xtknn=x]=(R1/n′f(x)4!-Rf(x))1n4,k∈{0,…,n-1}.
We clearly have ∥R1/n′f∥∞⩽CC8(V)l∥f∥8,∞. This, together with (6.4) and (6.1), completes the proof.
∎
Proof of Theorem 3.2.
We divide the proof into two steps.
Step 1. Let us prove that (2.4) is satisfied. We have
12(V12V0-V0V12)f(x)=(∂xV1(x)(V1(x)∂xV0(x)-∂xV1(x)V0(x))
+V1(x)(∂x2V0(x)V1(x)-∂x2V1(x)V0(x)))∂xf(x)
+V1(x)(V1(x)∂xV0(x)-∂xV1(x)V0(x))∂x2f(x).
Since V1(x)≠0, if we take
V~(x)=V1(x)(V1(x)∂xV0(x)-∂xV1(x)V0(x)),
then, using (3.3), V~ is well defined and satisfies (2.4).
Step 2. Now we are going to show the convergence in total variation distance. In order to do it we will use a result from [3]. First, applying the same reasoning as in the proof of Theorem 3.1, we can show that there exists some universal constants C,l⩾1 such that
(6.7)|〈g,Ptkn,tk+1nnf-Qtkn,tk+1nf〉|⩽n-4CC8(V)l∥g∥1,8∥f∥∞
with 〈⋅,⋅〉 the scalar product in L2(ℝ). Now we have (6.6) and (6.7), the result will be a consequence of
[3, Theorem 3.3], as soon as we check that the following ellipticity assumption holds:
there exists λ∗>0 such that infk⩽ninfx∈ℝ(∂w1ψk(x,w1,w0)|w1=w0=0)2⩾λ∗.
We fix k and we look at ψk(x,w1,w0) defined in (2.9). We suppose that ρk=3,ϵk=1 (the proof for ρk=1,2 or ϵk=-1 is similar). We consider the process xt(w~),0⩽t⩽T3, with Ti=i, w~=(w1,w0), solution of the following equation:
xt(w~)=x+w0∫0tV~(xs(w~))𝑑s,T0⩽t⩽T1,
xt(w~)=xT1(w~)+w1∫T1tV1(xs(w~))𝑑s,T1⩽t⩽T2,
xt(w~)=xT2(w~)+w0∫T2tV0(xs(w~))𝑑s,T2⩽t⩽T3.
We notice that ψk(x,wk+11,wk+10)=xT3(w~k+1) and consequently we have
∂zψk(x,wk+11,wk+10)=∂w1xT3(w~k+1).
Moreover, ∂w1xt(w)=0 for t⩽T1.
Now, let T1⩽t⩽T2. Then ∂w1xt(w~) solves the equation
∂w1xt(w~)=∫T1tV1(xs(w~))𝑑s+w1∫T1t∂V1(xs(w~))∂w1xs(w~)𝑑s.
It follows that
∂w1xt(w~)|w~=0=∫T1tV1(xs(0))ds=V1(x)(t-T1).
Notice that T2-T1=1. Then we have
∂w1xT3(w~)|w1=0=∂w1xT2(w~)|w~=0=V1(x).
and then, by (3.5),
(∂w1xT3(0))2⩾λ∗.
∎
References
[1]
Alfonsi A.,
High order discretization schemes for the CIR process: Application to affine term structure and Heston models,
Math. Comp. 79 (2010), no. 269, 209–237.
10.1090/S0025-5718-09-02252-2Search in Google Scholar
[2]
Bally V. and Clément E.,
Integration by parts formula and applications to equations with jumps,
Probab. Theory Related Fields 151 (2011), no. 3–4, 613–657.
10.1007/s00440-010-0310-ySearch in Google Scholar
[3]
Bally V. and Rey C.,
Approximation of Markov semigroup in total variation disctance,
preprint 2015, https://hal.archives-ouvertes.fr/hal.
10.1214/16-EJP4079Search in Google Scholar
[4]
Bally V. and Talay D.,
The law of the Euler scheme for stochastic differential equations. I: Convergence rate of the distribution function,
Probab. Theory Related Fields 104 (1996), no. 1, 43–60.
10.1007/BF01303802Search in Google Scholar
[5]
Bally V. and Talay D.,
The law of the Euler scheme for stochastic differential equations. II: Convergence rate of the density,
Monte Carlo Methods Appl. 2 (1996), no. 2, 93–128.
10.1515/mcma.1996.2.2.93Search in Google Scholar
[6]
Bossy M., Gobet E. and Talay D.,
A symmetrized Euler scheme for an efficient approximation of reflected diffusions,
J. Appl. Probab. 41 (2004), no. 3, 877–889.
10.1239/jap/1091543431Search in Google Scholar
[7]
Gobet E.,
Weak approximation of killed diffusion using Euler schemes,
Stochastic Process. Appl. 87 (2000), no. 2, 167–197.
10.1016/S0304-4149(99)00109-XSearch in Google Scholar
[8]
Gobet E. and Menozzi S.,
Stopped diffusion processes: Boundary corrections and overshoot,
Stochastic Process. Appl. 120 (2010), no. 2, 130–162.
10.1016/j.spa.2009.09.014Search in Google Scholar
[9]
Guyon J.,
Euler scheme and tempered distributions,
Stochastic Process. Appl. 116 (2006), no. 6, 877–904.
10.1016/j.spa.2005.11.011Search in Google Scholar
[10]
Jacod J., Kurtz T. G., Méléard S. and Protter P.,
The approximate Euler method for Lévy driven stochastic differential equations,
Ann. Inst. Henri Poincaré Probab. Statist. 41 (2005), no. 3, 523–558.
10.1016/j.anihpb.2004.01.007Search in Google Scholar
[11]
Jourdain B. and Kohatsu-Higa A.,
A Review of Recent Results on Approximation of Solutions of Stochastic Differential Equations,
Progr. Probab. 65,
Birkhäuser, Basel, 2011.
10.1007/978-3-0348-0097-6_9Search in Google Scholar
[12]
Kloeden P. E. and Platen E.,
Numerical Solution of Stochastic Differential Equations,
Appl. Math. (New York) 23,
Springer, Berlin, 1992.
10.1007/978-3-662-12616-5Search in Google Scholar
[13]
Kohatsu-Higa A. and Tankov P.,
Jump-adapted discretization schemes for Lévy-driven SDEs,
Stochastic Process. Appl. 120 (2010), no. 11, 2258–2285.
10.1016/j.spa.2010.07.001Search in Google Scholar
[14]
Konakov V. and Menozzi S.,
Weak error for stable driven stochastic differential equations: Expansion of the densities,
J. Theoret. Probab. 24 (2011), no. 2, 454–478.
10.1007/s10959-010-0291-xSearch in Google Scholar
[15]
Konakov V., Menozzi S. and Molchanov S.,
Explicit parametrix and local limit theorems for some degenerate diffusion processes,
Ann. Inst. Henri Poincaré Probab. Stat. 46 (2010), no. 4, 908–923.
10.1214/09-AIHP207Search in Google Scholar
[16]
Kusuoka S.,
Approximation of expectation of diffusion processes based on Lie algebra and Malliavin calculus,
Advances in Mathematical Economics. Vol. 6,
Adv. Math. Econ. 6,
Springer, Tokyo (2004), 69–83.
10.1007/978-4-431-68450-3_4Search in Google Scholar
[17]
Kusuoka S.,
Gaussian K-scheme: Justification for KLNV method,
Advances in Mathematical Economics. Vol. 17,
Adv. Math. Econ. 17,
Springer, Tokyo (2013), 71–120.
10.1007/978-4-431-54324-4_3Search in Google Scholar
[18]
Lyons T. and Victoir N.,
Cubature on Wiener space,
Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 460 (2004), no. 2041, 169–198.
10.1098/rspa.2003.1239Search in Google Scholar
[19]
Milstein G.,
Weak approximation of solutions of systems of stochastic differential equations,
Numerical Integration of Stochastic Differential Equations,
Math. Appl. 313,
Springer, Dordrecht (1995), 101–134.
10.1007/978-94-015-8455-5_4Search in Google Scholar
[20]
Ninomiya S. and Victoir N.,
Weak approximation of stochastic differential equations and application to derivative pricing,
Appl. Math. Finance 15 (2008), no. 1–2, 107–121.
10.1080/13504860701413958Search in Google Scholar
[21]
Protter P. and Talay D.,
The Euler scheme for Lévy driven stochastic differential equations,
Ann. Probab. 25 (1997), no. 1, 393–423.
10.1214/aop/1024404293Search in Google Scholar
[22]
Talay D. and Tubaro L.,
Expansion of the global error for numerical schemes solving stochastic differential equations,
Stochastic Anal. Appl. 8 (1990), no. 4, 483–509.
10.1080/07362999008809220Search in Google Scholar
![Figure 1 Log-Log representation of |𝔼[f(Xt)]-𝔼[f(Xtn)]|${|\mathbb{E}[f(X_{t})]-\mathbb{E}[f(X^{n}_{t})]|}$ for x=0.8${x=0.8}$, T=1${T=1}$, a=0.2${a=0.2}$, σ=2${\sigma=2}$, with respect to n for f(x)=𝟙x>1.1${f(x)=\mathbb{1}_{x>1.1}}$.](/document/doi/10.1515/mcma-2016-0120/asset/graphic/j_mcma-2016-0120_fig_x1.png)
