TY - JOUR
AU - Yamada, Toshihiro
AB - Abstract The paper proposes a new second-order weak approximation scheme for hypoelliptic diffusions or degenerate systems of stochastic differential equations satisfying a certain Hörmander condition. The scheme is constructed by a Gaussian process and a stochastic polynomial weight through a technique based on Malliavin calculus, and is implemented by a Monte Carlo method and a quasi-Monte Carlo method. A variance analysis for the Monte Carlo method is discussed, and further control variate methods are introduced to reduce the variance. The effectiveness of the proposed scheme is illustrated through numerical experiments for some hypoelliptic diffusions. 1. Introduction Hypoelliptic diffusions or degenerate systems of stochastic differential equations (SDEs) naturally appear in mathematical modelling in natural and social sciences such as physics, neural cell biology and financial mathematics. In particular, we typically encounter numerical problems of estimating probability distributions or computing expectations of hypoelliptic diffusions with nonsmooth test functions. On the other hand, constructing an efficient approximation scheme for these problems is not an easy task because we must pay attention to specific mathematical properties of hypoelliptic diffusions, especially degeneracy, see Kusuoka & Stroock (1985); Watanabe (1987); Ikeda & Watanabe (1989); Kloeden & Platen (1999); Milstein & Tretyakov (2004) and Hairer et al. (2011) for example. Following the study of the weak Euler–Maruyama scheme for SDEs by Bally & Talay (1996) under Hörmander’s condition and using a nonsmooth test function Kusuoka (2001, 2004) and Lyons & Victoir (2004) proposed a general methodology for higher order weak approximation called KLV (Kusuoka-Lyons-Victoir) method. Some algorithms based on the KLV method have been derived in the literature, for example, see Ninomiya & Victoir (2008); Gyurkó & Lyons (2010a,b); Litterer & Lyons (2012) and Bayer et al. (2013). In this paper we introduce a new higher order weak approximation scheme for |$N$|-dimensional SDE $$\begin{align} \textrm{d}X^x_t= {V}_0(X^x_t)\textrm{d}t+ \sum_{j=1}^d V_j(X^x_t) \textrm{d}W_t^{\,j}, \ \ \ X^x_0=x \in \mathbf{R}^{N}, \end{align}$$(1.1) driven by |$d$|-dimensional Brownian motion whose vector fields satisfy a type of Hörmander condition that naturally appears in mathematical modelling. For the model we introduce a second-order weak approximation using techniques from Malliavin calculus. The proposed scheme is inspired by a hypoelliptic example of Kolmogorov (1934). Let us consider the system: $$\begin{eqnarray} \mathrm{d} X_t^x= \sigma \textrm{d}W_t^1, \ \ \ \ \textrm{d}X^{x,2}_t= X_t^{x,1}\textrm{d}t, \end{eqnarray}$$(1.2) that is, an example of (1.1) with |$N=2$|⁠, |$d=1$| and |$V_0^1(x)=0$|⁠, |$V_1^1(x)=\sigma>0$|⁠, |$V_0^2(x)=x_1$|⁠, |$V_1^2(x)=0$| for |$x \in \mathbf{R}^{2}$|⁠. The system (1.2) is degenerate since |$V_1^2(x)=0$|⁠, but the solution of the SDE still has a density because Hörmander’s condition holds, i.e. |$V_1, [V_0,V_1]$| linearly span |$\mathbf{R}^2$| for any points on |$\mathbf{R}^2$|⁠. Our idea for the scheme is as follows: if we construct a sharp small time expansion of the expectation |$E[\,f(X_t^x) ]$| for the SDE (1.1) around the model $$\begin{eqnarray} \mathrm{d}\bar{X}^x_t= {V}_0(x)\textrm{d}t+ \sum_{j=1}^d V_j(x) \textrm{d}W_t^{\,j}, \ \ \ \bar{X}^x_0=x \in \mathbf{R}^{N}, \end{eqnarray}$$(1.3) which satisfies a certain Hörmander condition and has a density as in (1.2) we will have a higher order weak approximation for (1.1). Indeed, we show the following discretization for the system (1.1): $$\begin{eqnarray} E[\,f(X_T^x) ]= (Q_{T/n})^n f (x) + O \Big( \frac{1}{n^{2}} \Big), \end{eqnarray}$$(1.4) which will hold even if the test function |$f: \mathbf{R}^{N} \rightarrow \mathbf{R}$| is nonsmooth. Here, the operators |$\{ Q_t \}_{t>0}$| satisfy |$E[\,f(X_t^x) ]= Q_{t}\ f(x) + O (t^3)$| for small |$t$|⁠. In particular, the local approximation operator |$Q_t$| is given in the form |$(Q_{t}\ f)(x)=E[\,f(\bar{X}_t^x) \pi _t^x ]$| with a Malliavin weight|$\pi _t^x$| obtained by a certain polynomial of Brownian motions using integration by parts on Wiener space, by applying a perturbation expansion approach with Malliavin calculus. The advantages of using a perturbation expansion are as follows. First, the expansion based on Watanabe (1987) provides a refined, tailor-made small-time approximation around a certain nondegenerate Wiener functional, especially around the Kolmogorov example in this paper. Secondly, this approach enables us to avoid the simulation of the Lévy area term, which naturally appears in higher order schemes. In particular, through Malliavin integration by parts, we are able to get an explicit polynomial of Brownian motions for the Lévy area term as a Malliavin weight, which is an important feature from the viewpoint of implementation. Thirdly, the expansion may be regarded as a kind of extension of that of Bally & Talay (1996) in the sense that we use the Gaussian term (the Euler–Maruyama term) as the dominant part, and the Malliavin weight works as the correction term that contributes to higher order discretization. Finally, the expansion approach makes the treatment of nonsmoothness of the test functions easier and more flexible in weak approximations. Then, we can deal with many applications. Due to these reasons, the approach of the paper gives a powerful method and we will see its advantages in the numerical examples. In the error estimate for the discretization Malliavin’s integration by parts with Kusuoka–Stroock’s moment estimate for the inverse of the Malliavin covariance matrix (Kusuoka & Stroock, 1985) will play an important role. Then, the upper bound of the discretization error does not depend on any derivatives of the test function |$f$|⁠. This is the reason why the scheme works under nonsmooth test functions. The weak approximation is an extension of Takahashi & Yamada (2016); Yamada & Yamamoto (2018a,b); Yamada (2019) in the sense that the model is extended to the hypoelliptic case. We perform numerical experiments in order to check the validity of the proposed scheme using some important hypoelliptic models from applied mathematics. In particular, we give an Asian option approximation using a second-order scheme that is an extension of the methods of Kunitomo & Takahashi (1992) and Takahashi (2015), and show comparison with the method of Ninomiya & Victoir (2008). Furthermore, we give a nontrivial multidimensional example for pricing volume-weighted average price (VWAP) option, which is an interesting recent example arising in financial mathematics. To the best of our knowledge it is the first attempt at deriving a higher order discretization scheme for pricing a VWAP option. The organization of the paper is as follows. In Section 2, after we introduce basic notations and mathematical conditions, we discuss the construction of our second-order weak approximation scheme (the main result is given as Theorem 2.3). Section 3 shows the order of the variance of the proposed scheme from a mathematical aspect. In Section 4 we provide numerical examples to show the validity of our new scheme and, in particular, we compare the new second-order scheme to another second-order scheme in the literature. Section 5 concludes the paper and some technical lemmas are proved in the Appendix. 2. Weak approximation We introduce some notations. Let |${C}_b^{\infty }(\mathbf{R}^{m},\mathbf{R}^{n})$| be a space of infinitely differentiable functions |$f: \mathbf{R}^m \to \mathbf{R}^n$| such that |$f$| and its partial derivatives are bounded. We write |${C}_b^{\infty }(\mathbf{R}^{m})$| for |${C}_b^{\infty }(\mathbf{R}^{m},\mathbf{R})$|⁠. For a bounded and measurable function |$f: \mathbf{R}^{m} \rightarrow \mathbf{R}$| we define |${\|\,f \|_{\infty }:= \sup _{x \in \mathbf{R}^{m}}|f(x)|}$|⁠, for |$f \in{C}_b^{\infty }(\mathbf{R}^{m})$|⁠, and define $$\begin{align*} \|\nabla^{k} f \|_{\infty}:= \max_{j_1,\ldots, j_k\in \{1,\ldots,N\}} \bigg\| \frac{\partial^{k}f}{\partial x_{j_1}\cdots \partial x_{j_k}} \bigg\|_{\infty}. \end{align*}$$ We may write |$\partial _i \varphi (x)$| in short instead of |${\frac{\partial \varphi }{ \partial x_i}(x)} $| for |$\varphi \in{C}_b^{\infty }( \mathbf{R}^m)$|⁠. For all |$z \in \mathbf{R}^m$| we write |${|z|: = \sqrt{ \sum _{i=1}^m |z_i|^2}}$|⁠. 2.1 Malliavin calculus We introduce some basic results and notations of Malliavin calculus. See Ikeda & Watanabe (1989); Shigekawa (2004) and Nualart (2006) for more details. Let |$( \varOmega , \mathcal{F}, \mathbb{P})$| be |$d$|-dimensional Wiener space, where |$\varOmega =\{ w: [0,T] \to \mathbf{R}^d ; w \;\textrm{is continuous and}\;w_0=0 \}$| is a Banach space induced by the uniform norm on |$[0,T]$|⁠, |$\mathcal{F}$| is the Borel sigma-field over |$\varOmega $| and |$\mathbb{P}$| is the Wiener measure such that the canonical process |$\{ W_t \}_{t \in [0,T]}$| is a |$d$|-dimensional Brownian motion. We will denote by |$H$| the Hilbert space |$L^2 ( [0,T], \mathbf{R}^d )$| with the inner product |$\langle \cdot , \cdot \rangle _H$| given by |${\langle h_1,h_2 \rangle _H=\int _0^T h_1(s) \cdot h_2(s)ds}$| for |$h_1,h_2 \in H$| where |$\cdot $| represents the Euclidean inner product. For |$h \in H$| we will use the notation |$W(h)={\int _0^T} h(s) \cdot \textrm{d}W_s$| for the Wiener integral. We write |$L^p (\varOmega , \mathbf{R}^N), 1 \leqslant p < \infty $| for a Banach space of measurable functions |$F: \varOmega \to \mathbf{R}^N$| such that |$\| F \|_p:= E [ |F|^p ]^{1/p} < \infty $|⁠. Let |$\mathcal{S}$| be the set of cylindrical random variables |$\mathcal{S}=\{ F= f ( W(h_1), \cdots , W(h_n) ); \ n \geqslant 1, \ f \in{C}_b^{\infty }( \mathbf{R}^n), \ h_1, \ldots , h_n \in H\}$|⁠. For |$F \in \mathcal{S}$| the derivative of |$F$| is the |$\mathbf{R}^d$|-valued stochastic process |$\{ D_t F\}_{t \in [0,T]}$| given by $$\begin{align} D_t F = \sum_{i=1}^n \partial_i f ( W(h_1), \cdots, W(h_n) ) h_i(t), \end{align}$$(2.1) and for |$j=1, \ldots , d$| the |$j$|th element of |$D_t F$| will be denoted by |$D_t^{\,j} F$|⁠. The |$k$|th order derivative can be also defined as a random vector on |$[0, T]^k \times \varOmega $| whose coordinates are given by |$D_{t_1,\cdots , t_k}^{j_1, \ldots , j_k} F = D_{t_k}^{j_k} \cdots D_{t_1}^{j_1} F$| for |$k \in \mathbf{N}$|⁠. One can show that the operator |$D^k$| is closable from |$\mathcal{S}$| into |$L^p (\varOmega , H^{\otimes k} )$| for any |$p \geqslant 1$|⁠. Then, for any |$k \in \mathbf{N}$| and |$p \geqslant 1$|⁠, we denote by |$\mathbf{D}^{k,p}$| the completion of |$\mathcal{S}$| with respect to the norm $$\begin{align} \| F \|_{k,p} = \Big\{ E [|F |^p] + \sum{}_{i=1}^k E[ \| D^i F \|_{H^{\otimes i}} ^p ] \Big\}^{1/p}. \end{align}$$(2.2) Furthermore, we define |$\mathbf{D}^{\infty } = {\bigcap _{k \in \mathbf{N},p \geqslant 1}} \mathbf{D}^{k,p}$| and its dual |$ {\mathbf{D}}^{- \infty } = (\mathbf{D}^{\infty })^{\prime}$|⁠. We denote by |${}_{-\infty } \langle \varPhi ,G \rangle _{{\infty }}$| the coupling between |$\varPhi \in \mathbf{D}^{-\infty }$| and |$G \in{\mathbf D}^{\infty }$|⁠. The adjoint operator of |$D$| is defined by |$\delta : L^2 ( [0, T] \times \varOmega ) \to L^2 ( \varOmega )$| through the duality formula $$\begin{align} E[ F \delta(u) ] = E[ \langle DF, u \rangle_H] \end{align}$$(2.3) for |$F \in \mathbf{D}^{1,2}$|⁠. In particular, if the stochastic process |$u$| is square integrable and is adapted to the filtration generated by |$\{W_t\}_{t \in [0,T]}$|⁠, |$u$| belongs to the domain of |$\delta $| and |$\delta (u)$| coincides with the Itô integral, i.e. the Skorokhod integral |$\delta (u) = {\int _0^T u_t \cdot \textrm{d}W_t}$|⁠. For a random vector |$F = (F^1, \ldots , F^m) \in ( \mathbf{D}^{\infty })^m$| the Malliavin covariance is defined by a symmetric non-negative definite |$m \times m$| matrix |$\sigma ^F_{i\,j} = \langle DF^i, DF^{\,j} \rangle _H, \ i,j =1, \ldots ,m.$| In particular, |$F$| is said to be nondegenerate in the Malliavin sense (hereafter, we simply say nondegenerate) if |$\sigma ^F$| is invertible a.s. and |$ (\det \sigma ^F)^{-1} \in {\bigcap _{p \geqslant 1}}L^p (\varOmega )$|⁠. The following is the integration by parts on the Wiener space, which often appears in this paper. Proposition 2.1 (Integration by parts). Let |$F \in (\mathbf{D}^{\infty })^m $| be a nondegenerate random vector. Denote the inverse of the Malliavin covariance by |$\gamma ^F =(\gamma ^F_{i\,j})_{i,j=1,\ldots , m}$|⁠. Also let |$G \in{\mathbf D}^\infty $| and |$ \varphi \in{C}_b^{\infty }(\mathbf{R}^m)$| be given. Then, for any multi-index |$\alpha =(\alpha _1,\ldots ,\alpha _k) \in \{1,2,\ldots ,m\}^k, k \geqslant 1$|⁠, there exists |$H_{\alpha }(F,G)$| such that $$\begin{equation} {E}[\partial^{\alpha} \varphi (F) G] = {E} [\varphi(F) H_{\alpha}(F,G)], \nonumber \end{equation}$$ where |$H_{(\alpha _1,\ldots ,\alpha _k)}(F,G)$| is determined recursively by $$\begin{align} H_{(i)}(F,G) = \sum_{j=1}^d \delta \left( \gamma^F_{i\,j} G DF^{\,j} \right), \ i=1,\cdots,N, \ \ \ H_{(\alpha_1,\ldots,\alpha_k)}(F,G) =H_{(\alpha_k)}(F, H_{(\alpha_1,\ldots,\alpha_{k-1})}(F,G)). \nonumber \end{align}$$ Let |$\mathcal{S}^{\prime}(\mathbf{R}^m)$| be the space of Schwartz distributions that is the dual of Schwartz space of rapidly decreasing functions |$f:\mathbf{R}^m \rightarrow \mathbf{R}$|⁠. It holds that the composition |$T(F)$| of |$T \in \mathcal{S}^{\prime}(\mathbf{R}^m)$| and a nondegenerate |$F \in (\mathbf{D}^{\infty })^m$| is well-defined as an element of |$ {\mathbf D}^{-\infty }$| and we have $$\begin{align} {}_{ {-\infty}} \langle \partial^{\alpha} T (F), G\rangle_{{\infty}} = {}_{ {-\infty}} \langle T(F), H_{\alpha}(F,G) \rangle_{{\infty}} \end{align}$$(2.4) for any multi-index |$\alpha =(\alpha _1,\ldots ,\alpha _k) \in \{1,2,\ldots ,m\}^k, k \geqslant 1$| and |$G \in{\mathbf D}^{\infty }$|⁠. Let us denote the Dirac delta function mass at |$y \in \mathbf{R}^m$| by |$\delta _y \in \mathcal{S}^{\prime}(\mathbf{R}^m)$|⁠. Then, for a nondegenerate |$F \in (\mathbf{D}^{\infty })^m$|⁠, the composition |$\delta _y (F)$| is in |$ {\mathbf D}^{-\infty }$| and one has $$\begin{align} E [\varphi(F) G] = \int_{\mathbf{R}^m} \varphi(y) {}_{ {-\infty}} \langle \delta_y( F ), G \rangle_{{\infty}} \textrm{d}y, \end{align}$$(2.5) for any |$\varphi \in C_b^{\infty } (\mathbf{R}^m)$| and |$G \in{\mathbf D}^{\infty }$|⁠. See Takahashi & Yamada (2012) and Takahashi (2015) for the application of generalized Wiener functionals, for example. 2.2 Stochastic differential equation On the |$d$|-dimensional Wiener space we consider the following |$N$|-dimensional Stratonovich-type SDE driven by |$d$|-dimensional Brownian motion: $$\begin{equation} \textrm{d}X^x_t= \widetilde{V}_0(X^x_t)\textrm{d}t+ \sum_{j=1}^d V_j(X^x_t) \circ \textrm{d}W_t^{\,j}, \ \ \ X^x_0=x \in \mathbf{R}^{N}, \end{equation}$$(2.6) and the corresponding Itô-type SDE is given by $$\begin{align} \textrm{d}X^x_t= {V}_0(X^x_t)\textrm{d}t+ \sum_{j=1}^d V_j(X^x_t) \textrm{d}W_t^{\,j}, \ \ \ X^x_0=x \in \mathbf{R}^{N}. \end{align}$$(2.7) In particular, we suppose that each component of (2.7) is given as follows. $$\begin{align} \textrm{d}X_t^{x,i_1}=&\; {V}_0^{i_1}( X^x_t) \textrm{d}t + \sum_{j=1}^d V_j^{i_1} (X^x_t) \textrm{d}W_t^{\,j}, \ \ \ i_1 = 1,2, \ldots, K, \end{align}$$(2.8) $$\begin{align} \textrm{d}X_t^{x,i_2}=&\; {V}_0^{i_2}(X^x_t)\textrm{d}t, \ \ \ i_2 = K+1, \ldots, N, \end{align}$$(2.9) where |$K \in \{1, 2, \ldots , N\}$|⁠. We assume that |$\widetilde{V}_0, V_j$|⁠, |$j=1,\ldots ,d$| satisfy |$\widetilde{V}_0, V_j \in{C}^{\infty }_b (\mathbf{R}^N,\mathbf{R}^N)$|⁠, which are identified with vector fields through |$\widetilde{V}_0 ={\sum _{i=1}^{N} \widetilde{V}_0^i(x)} \frac{\partial }{\partial x_i}$| and |$V_j = {\sum _{i=1}^N V_j^i(x) \frac{\partial }{\partial x_i}}.$| For |$V, W \in{C}^{\infty }_b (\mathbf{R}^N,\mathbf{R}^N)$| the Lie bracket is defined by |$[V,W] = VW - WV$|⁠. Hereafter, we use two types of differential operators for stochastic Taylor expansion: for |$\varphi \in C_b^{\infty }(\mathbf{R}^N)$| and |$x \in \mathbf{R}^N$| $$\begin{align} \hat{V}_0 \varphi (x) =&\; \sum_{i=1}^N V_0^i (x)\frac{\partial \varphi}{\partial x_i}(x) + \frac{1}{2} \sum_{i,j=1}^N \sum_{k=1}^d V_k^i(x) V_k^{\,j}(x) \frac{\partial^2 \varphi}{\partial x_i \partial x_{\,j}}(x), \end{align}$$(2.10) $$\begin{align} \hat{V}_j \varphi (x) =&\; \sum_{i=1}^N V_j^i (x)\frac{\partial \varphi}{\partial x_i}(x), \ \ j=1, \ldots, d, \end{align}$$(2.11) and also for a multi-index |$\alpha = (\alpha _1, \ldots , \alpha _l)\in \{0, 1, \ldots , d \}^\ell $|⁠, |$\ell \in \mathbf{N}$|⁠, we define the iterated Itô integral as $$\begin{align*} I_{\alpha}(t) = \int_{0 < t_1 < \cdots < t_\ell < t} \textrm{d}W_{t_1}^{\alpha_1} \cdots \textrm{d}W_{t_\ell}^{\alpha_\ell}, \ \ \ \ t \geqslant 0. \end{align*}$$ For a multi-index |$\alpha = (\alpha _1, \ldots , \alpha _\ell )\in \{0, 1, \ldots , d \}^\ell $|⁠, |$\ell \in \mathbf{N}$|⁠, define |$|\alpha |=\ell $| and |$\| \alpha \|=\#\{ j ; \alpha _j\neq 0 \}+2 \#\{ j ; \alpha _j= 0 \}=| \alpha |+\#\{ j ; \alpha _j= 0 \}$|⁠. Furthermore, we introduce an |$(N-K)$|-dimensional vector |$\bar{V}_j (x) = ( \hat{V}_j V_0^{K+1}(x), \ldots , \hat{V}_j V_0^N(x) )$| for |$j=1,2, \ldots ,d$| and |$x \in \mathbf{R}^N$|⁠. Throughout the paper we put the following condition on the coefficients of the SDE. Assumption 2.2 |$V_1(x), \ldots , V_d(x)$| linearly span |$\mathbf{R}^K$| for all |$x \in \mathbf{R}^N$|⁠. |$\bar{V}_1(x), \ldots , \bar{V}_d(x)$| linearly span |$\mathbf{R}^{N-K}$| for all |$x \in \mathbf{R}^N$|⁠. Noting that |$[V_j, \widetilde{V}_0 ]^i (x) = \hat{V}_j V_0^i(x)$| for |$i=K+1, \ldots ,N$| and |$j=1, \ldots ,d$| we see that under Assumption 2.2, $$\begin{align} V_1(x), V_2(x), \ldots, V_d(x), [V_1, \widetilde{V}_0](x), [V_2, \widetilde{V}_0](x),\ldots, [V_d, \widetilde{V}_0](x) \end{align}$$(2.12) linearly span |$\mathbf{R}^N$| for every |$x \in \mathbf{R}^N$|⁠. We denote by |$\{P_t\}_{t \geqslant 0}$| the semigroup of linear operators defined by $$\begin{equation} (P_t\ f)(x)={E} [\,f(X_t^x)], \ \ \ t\geqslant 0, \ x \in \mathbf{R}^{N}, \end{equation}$$(2.13) for a bounded measurable function |$f:\mathbf{R}^{N} \rightarrow \mathbf{R}$|⁠. Our goal in the paper is to construct a second-order discretization of |$(P_Tf)(x)$| for |$T\geqslant 1$|⁠, |$x \in \mathbf{R}^{N}$| and a bounded measurable function |$f:\mathbf{R}^{N} \rightarrow \mathbf{R}$|⁠. Here, we summarize the strategy of the construction. Step 1. Local approximation: Using the perturbation expansion method we obtain local approximation operators |$\{Q_t\}_{t>0}$| satisfying $$\begin{align} P_t \varphi(x) = Q_t \varphi (x) + O(t^3) \end{align}$$(2.14) for |$t \in (0,1], \; x \in \mathbf{R}^N$| and |$\varphi \in C_b^{\infty }(\mathbf{R}^N)$|⁠. We will discuss the local approximation in Theorem 2.3 in Section 2.3 and describe the derivation of |$Q_t$| in Subsection 2.3.1. In particular, the error |$O(t^3)$| is decomposed into two terms |$\mathcal{E}_1^{\varphi } (t,x)$| and |$\mathcal{E}_2^{\varphi } (t,x)$|⁠, which are obtained by Lemmas 2.4 and 2.6. Step 2. Global approximation (main result): Let |$T \geqslant 1$|⁠. We compose the local approximation operator |$n$|-times as $$\begin{eqnarray*} (Q_{T/n})^n f (x) = (\underbrace{Q_{T/n} \circ \cdots \circ Q_{T/n} \circ Q_{T/n}}_{n{\text{-times}}} f )(x), \ \ \ x \in \mathbf{R}^N, \end{eqnarray*}$$ and show that $$\begin{align} P_T f (x) = (Q_{T/n})^n f (x) + O (n^{-2}), \end{align}$$(2.15) where |$f:\mathbf{R}^N \to \mathbf{R}$| is bounded and measurable (not necessarily smooth). This is given as our main result in Section 2.4, and its proof is shown in Section 2.4.1 with technical lemmas some of which are proved in the Appendix. 2.3 Local approximation First, we give a local approximation for |$(P_t\varphi )(x)$|⁠, |$\varphi \in C_b^{\infty } (\mathbf{R}^N)$|⁠, of the order |$O(t^3)$| to obtain the discretization for |$(P_Tf)(x)$| of the order |$O(n^{-2})$|⁠. We apply a perturbation expansion approach with Malliavin calculus to obtain a refined, tailor-made local approximation, which can be explicitly implemented by the Monte Carlo method. Let |$ \bar{X}_t^{x}=( \bar{X}_t^{x,1},\cdots , \bar{X}_t^{x,N})$|⁠, |$t\geqslant 0$| be a process given by $$\begin{eqnarray} \bar{X}_t^{x,i_1}&=&x_{i_1} + {V}_0^{i_1}(x)t + \sum_{j=1}^d V_j^{i_1}(x) W^{\,j}_t, \ \ i_1=1,\ldots,K, \end{eqnarray}$$(2.16) $$\begin{eqnarray} \bar{X}_t^{x,i_2}&=& x_{i_2} + {V}_0^{i_2}(x)t + \sum_{j=0}^{d} \hat{V}_j {V}_0^{i_2}(x) \int_{0}^{t}W^{\,j}_{s} \textrm{d}s, \ \ i_2=K+1, \ldots, N, \end{eqnarray}$$(2.17) where |$W_t^0 = t$|⁠. Under Assumption 1 it is guaranteed that for every |$i_2$| there exists at least one |$j \in \{1, 2,\ldots ,d\}$| such that |$\hat{V}_{j}V_0^{i_2}(x) \neq 0$|⁠, namely, each |$\bar{X}_t^{i_2}(x)$| has stochastic term |$ \hat{V}_j V_0^{i_2}(x) {\int _0^t W_s^{\,j} \textrm{d}s}$|⁠. Then, we introduce a set |${I:= \bigcup _{m=K+1}^N} \{ \alpha \in \{ 1, 2,\ldots , d \} \: \ [V_{\alpha }, V_0 ]^m (\cdot ) \neq 0 \}$|⁠, where |$[V_{\alpha }, V_0 ]^m (x)$| means the |$m$|th component of |$N$|-dimensional vector |$[V_{\alpha }, V_0 ](x)$| for |$x \in \mathbf{R}^N$|⁠. Moreover, for arbitrary fixed |$t>0$|⁠, |$\bar{X}_t^{x}$| also satisfies Hörmander condition and has a density under Assumption 1 (as in the example (1.2) of Kolmogorov, 1934). The next theorem gives the local approximation for |$P_t\ f(x)$| satisfying |$P_t\ f(x)=Q_{t}\ f(x)+O(t^3)$|⁠. Theorem 2.3 (Local approximation). There exist |$C>0$| and |$m \in \mathbf{N}$| such that $$\begin{align*} \sup_{x \in \mathbf{R}^N} \Big| P_t \varphi(x)- Q_t \varphi(x) \Big| \leqslant C \left( \| \varphi \|_{\infty} + \sum_{i=1}^m \| \nabla^i \varphi \|_{\infty} \right) t^3 \end{align*}$$ for all |$\varphi \in C_b^{\infty } (\mathbf{R}^N)$| and |$t>0$|⁠, where |$\{ Q_{t} \}_{t>0}$| are operators given by $$\begin{equation*}(Q_{t} \varphi)(x)= {E} [ \varphi (\bar{X}_t^x ) {\pi}_t^x ], \ \ \ (t,x) \in (0,T] \times \mathbf{R}^{N}, \ \varphi \in C_b^{\infty}(\mathbf{R}^{N}),\end{equation*}$$ with a weight $$\begin{eqnarray} {\pi}_t^x &=& 1 + \sum_{ i=1}^{K} \sum_{0 \leqslant j_1,j_2 \leqslant d} \hat{V}_{j_1} V_{j_2}^i(x) H_{(i)} \left( \bar{X}_t^x, \frac{1}{2}{W_t^{j_1}W_t^{j_2}}- \frac{t}{2} \mathbf{1}_{ j_1 = j_2 \neq 0} \right) \end{eqnarray}$$(2.18) $$\begin{eqnarray} &&+ \sum_{1 \leqslant i_1,i_2 \leqslant K} \sum_{1 \leqslant j_1,j_2 \leqslant d} \hat{V}_{j_1} V_{j_2}^{i_1}(x) \hat{V}_{j_1} V_{j_2}^{i_2}(x) H_{(i_1,i_2)}( \bar{X}_t^x, 1) \frac{t^2}{4} \end{eqnarray}$$(2.19) $$\begin{eqnarray} &&+ \sum_{i=1}^{K} \sum_{\alpha \in I} \sum_{\substack{0 \leqslant j \leqslant d\\ j \neq \alpha} } [ V_{\alpha}, V_j ]^i(x) H_{(i)} \left( \bar{X}_t^x, \frac{1}{t} \left(\int_0^t W_s^{\alpha} \textrm{d}s \right) W_t^{j}- \frac{1}{2} W_t^{\alpha} W_t^{\,j} \right). \end{eqnarray}$$(2.20) In particular, we have the error decomposition $$\begin{eqnarray*} P_t \varphi(x)- Q_t \varphi(x)= \mathcal{E}_1^{\varphi} (t,x)+ \mathcal{E}_2^{\varphi} (t,x), \ \ \ (t,x) \in (0,T] \times \mathbf{R}^{N},\end{eqnarray*}$$ where the error terms satisfy: there is |$C>0$| such that $$\begin{equation*}\sup_{x \in \mathbf{R}^N} | \mathcal{E}_1^{\varphi}(t,x) | \leqslant C \| \varphi \|_{\infty}t^3 \end{equation*}$$ for |$\varphi \in C_b^{\infty } (\mathbf{R}^N)$| and |$t>0$|⁠, and $$\begin{equation*}\mathcal{E}_2^{\varphi}(t,x)=t^3 \sum_{\lambda_l,l\leqslant m} E[ \partial^{\lambda_l} \varphi(\bar{X}_t^x)] {\mathcal G}_{\lambda_l}(t,x),\end{equation*}$$ for some |$m\geqslant 1$|⁠, |$\lambda _j \in \{1,\cdots ,N \}^{\,j}$|⁠, |$j\geqslant 1$|⁠, and function |${\cal G}_{\lambda _l}$| such that |${ \cal G}_{\lambda _l}(\cdot ,x)$| is nondecreasing for any fixed |$x \in \mathbf{R}^N$| and |${\cal G}_{\lambda _l}(t,\cdot)$| is bounded for any fixed |$t \in (0,1]$|⁠. 2.3.1 Proof of Theorem 2.3 Consider the following perturbed SDE: $$\begin{eqnarray*} \textrm{d}X_t^{x,i_1,\varepsilon}&=&\varepsilon^2 V_0^{i_1}(X_t^{x,\varepsilon})\textrm{d}t+\varepsilon \sum_{j=1}^d V_j^{i_1}(X_t^{x,\varepsilon})\textrm{d}W^{j}_t, \ \ \ i_1=1,\cdots,K, \\ \textrm{d}X_t^{x,i_2,\varepsilon}&=&\varepsilon^2 V_0^{i_2}(X_t^{x,\varepsilon})\textrm{d}t, \ \ \ i_2=K+1,\cdots,N, \end{eqnarray*}$$ where |$\varepsilon \in (0,1]$|⁠. Note that |$X_1^{x,\varepsilon }|_{\varepsilon =\sqrt{t}} \overset{\mathrm{Law}}{\sim } X_t^{x,1}=X_t^{x}$| in probability law. Then it is enough to show the expansion of |$E[\varphi (X_1^{x,\varepsilon })]$| of the order |$O(\varepsilon ^6)$| to get the local approximation for |$(P_tf)(x)$| of the order |$O(t^3)$|⁠. Let |$Y_t^{0,i_1}=\sum _{j=1}^d V_j^{i_1}(x)W^{j}_t$|⁠, |$Y_t^{0,i_2}=\sum _{j=1}^d \hat{V}_{j}V_0^{i_2}(x)\int _0^t W^{j}_s \textrm{d}s$|⁠, and for |$m\in \mathbf{N}$|⁠, we also introduce $$\begin{eqnarray} Y_t^{m,i_1}&=&\sum_{J=(j_1,\cdots,j_l),\|J\|=m+1,|J|\geqslant 2} \hat{V}_{j_{1}} \cdots \hat{V}_{j_{l-1}}V_{j_l}^{i_1}(x)I_{(j_1,\cdots,j_l)}(t), \end{eqnarray}$$(2.21) $$\begin{eqnarray} Y_t^{m,i_2}&=&\sum_{J=(j_1,\cdots,j_l), \|J\|=m+1,|J|\geqslant2} \hat{V}_{j_{1}} \cdots \hat{V}_{j_{l}}V_{0}^{i_1}(x) I_{(j_1,\cdots,j_l,0)}(t). \end{eqnarray}$$(2.22) Let |$\tau : \mathbf{R}^N \to \mathbf{R}^N$| be a map given by $$\begin{eqnarray} \tau^{i_1} (y) = &\;x_{i_1} + \varepsilon^2 V_0^{i_1}(x) + \varepsilon y_{i_1}, \; i_1 = 1, \ldots, K , \end{eqnarray}$$(2.23) $$\begin{eqnarray} \tau^{i_2}(y) =&\; x_{i_2} + \varepsilon^2 V_0^{i_2}(x) + \frac{\varepsilon^4}{2} \hat{V}_0 V_0^{i_2}(x) + \varepsilon^3 y_{i_2}, \; i_2 =2, \ldots, K+1, \ldots, N . \end{eqnarray}$$(2.24) We define a random vector |$\bar{X}_1^{x, \varepsilon } =( \bar{X}_1^{x,1, \varepsilon }, \ldots , \bar{X}_1^{x,N, \varepsilon })$| by |$\bar{X}_1^{x, i, \varepsilon } = \tau ^i(Y_1^0)$| for |$i=1,2, \ldots , N$|⁠. The following expansion is fundamental for our discussion. Lemma 2.4 There exists a |$C>0$| such that $$\begin{eqnarray*} &&\sup_{x \in \mathbf{R}^N} \Bigg| E[\varphi(X_1^{x,\varepsilon})]-\Bigg\{E[\varphi(\bar{X}_1^{x,\varepsilon})]\\ && +\sum_{j=1}^5 \varepsilon^{\,j} \sum_{k=1}^{\,j} \sum_{ \substack{ (\beta_i)_{i=1}^k \subset \mathbf{N} \\ s.t. \ \sum_{i=1}^k \beta_i=j}} \sum_{ (\alpha_1,\cdots,\alpha_k) \in \{1,\cdots,N \}^k } \frac{1}{k!} E\left[ \varphi(\bar{X}_1^{x,\varepsilon}) H_{(\alpha_1,\cdots,\alpha_k)}\left(Y_1^{0},\prod_{l=1}^k Y_1^{\beta_l,\alpha_l}\right) \right] \Bigg\} \Bigg|\\ &&\leqslant C \| \varphi \|_{\infty} \varepsilon^6, \end{eqnarray*}$$ for all |$\varphi \in C_b^{\infty } (\mathbf{R}^N)$| and |$\varepsilon \in (0,1]$|⁠. Proof of Lemma2.4. See Appendix A. Some expansion coefficients appearing in Lemma 2.4 will be also |$O(\varepsilon ^6)$|⁠. We note that the product |${\prod _{l=1}^k} Y_1^{\beta _l,\alpha _l}$| in Lemma 2.4 can be transformed into a sum of single iterated stochastic integrals using Itô’s formula. Indeed, for |$j=1,\cdots ,5$|⁠, |$k=1\cdots ,j$|⁠, |$\beta _i \geqslant 1$|⁠, |$i=1,\cdots ,k$|⁠, such that |$\sum _{i=1}^k \beta _i=j$| and a multi-index |$(\alpha _1,\cdots ,\alpha _k) \in \{1,\cdots ,N \}^k$|⁠, we have $$\begin{align*} \prod_{l=1}^k Y_1^{\beta_l,\alpha_l}=\sum_{r=1,\cdots,m} c_{r}(x) I_{\nu^r}(1) \end{align*}$$ for some |$m \in \mathbf{N}$|⁠, |$c_r(x) \in \mathbf{R}$| and multi-indices |$\nu ^r \in \{ 0,1,\cdots ,d \}^{e}$| such that |$\| \nu ^r \|=j+k$|⁠, |$r=1,\cdots ,m$|⁠. Then, the expansion terms in Lemma 2.4 can be written as $$\begin{align} \varepsilon^{\,j} E\left[ \varphi(\bar{X}_1^{x,\varepsilon}) H_{(\alpha_1,\cdots,\alpha_k)}\left(Y_1^{0},\prod_{l=1}^k Y_1^{\beta_l,\alpha_l}\right) \right]=&\; \varepsilon^{j+\iota( \alpha)} E\left[ \partial^{\alpha} \varphi(\bar{X}_1^{x,\varepsilon}) \prod_{l=1}^k Y_1^{\beta_l,\alpha_l} \right] \end{align}$$(2.25) $$\begin{align} = &\; \sum_{r=1,\cdots,m} c_{r}(x) \varepsilon^{\iota(\alpha)-k} E[ \partial^{\alpha} \varphi(\bar{X}_1^{x,\varepsilon}) \varepsilon^{j+k} I_{\nu^r}(1) ], \end{align}$$(2.26) where |$\iota (\alpha )$| is an integer depending on |$\alpha $| such that |$\iota ( \alpha ) \geqslant k$|⁠. In (2.25) we used the elementary integration by parts: $$\begin{align} E[ \varphi(\bar{X}_1^{x,\varepsilon}) H_{(i_1)}(Y_1^0,G) ]=&\;\varepsilon E[ \partial_{i_1} \varphi(\bar{X}_1^{x,\varepsilon}) G ], \ \ i_1=1,\cdots,K, \end{align}$$(2.27) $$\begin{align} E[ \varphi(\bar{X}_1^{x,\varepsilon}) H_{(i_2)}(Y_1^0,G) ]=&\; \varepsilon^3 E[ \partial_{i_2} \varphi(\bar{X}_1^{x,\varepsilon}) G ], \ \ i_2=K+1,\cdots,N, \end{align}$$(2.28) for |$G \in \mathbf{D}^{\infty }$|⁠. For the type of expectations in (2.26) we have the following result. Lemma 2.5 Let |$m \geqslant 3$| and |$\beta = (\beta _1, \ldots , \beta _m) \in \{0,1, \ldots , d\}^m$| be a multi-index. Define a multi-index |$\hat{\beta } = ( \beta _{i_1}, \ldots , \beta _{i_e}) \in \{ 1, 2, \ldots , d\}^e$| where |$0 \leqslant e \leqslant m$| and |$1 \leqslant i_1 < \cdots < i_e \leqslant m$| by collecting all nonzero elements in |$\beta $|⁠. Then, it holds that for |$\varphi \in C^{\infty }_b (\mathbf{R}^N)$| and |$\gamma \in \{1, \ldots , N\}^l$|⁠, |$l \in \mathbf{N}$|⁠, $$\begin{eqnarray} &&E \big[ \partial^{\gamma} \varphi(\bar{X}_1^{x,\varepsilon}) \varepsilon^{\| \beta \|} I_{\beta}(1) \big] \end{eqnarray}$$(2.29) $$\begin{eqnarray} &&= \sum_{\alpha \in \{1, 2, \ldots, N \}^e} \varepsilon^{\| \beta \|} E \big[ \partial^{\alpha} \partial^{\gamma} \varphi(\bar{X}_1^{x,\varepsilon}) \big] \int_{0 < t_{1} < \cdots < t_{m} < 1} {\psi}_{\beta_{i_1}}^{\alpha_1,\varepsilon}(t_{{i_1}},x) \cdots{\psi} _{\beta_{i_e}}^{\alpha_e,\varepsilon}(t_{i_e},x) \textrm{d}t_{1} \cdots \textrm{d}t_{m}, \end{eqnarray}$$(2.30) where the deterministic function |${\psi }_{j}^{i,\varepsilon }(s,x), \; (s,x) \in [0,T] \times \mathbf{R}^N$| for |$ i=1, \ldots , N,\ j=0, \ldots , d$|⁠, and |$\varepsilon \in (0,1]$| is given by $$\begin{equation*} {\psi}_{j}^{i,\varepsilon}(s,x) = \begin{cases} \varepsilon V_j^i(x) \textbf{1}_{0 \leqslant s \leqslant 1} & \textrm{if} \ \ i = 1, \ldots, K\\ \varepsilon^3 \hat{V}_j V_0^i(x) (1-s) \textbf{1}_{0 \leqslant s \leqslant 1} & \textrm{if} \ \ i = K+1, \ldots, N. \end{cases} \end{equation*}$$ Furthermore, there are an integer |$e \geqslant 1$| and a constant |$C> 0$| such that $$\begin{align} {\sup_{x \in \mathbf{R}^N}} \left| E[ \partial^{\gamma} \varphi( \bar{X}_1^{x,\varepsilon} ) \varepsilon^{\| \beta \|} I_{\beta}(1) ] \right| \leqslant \varepsilon^6 C \| \nabla^e \varphi \|_{\infty} \end{align}$$(2.31) for all |$\varphi \in C^{\infty }_b (\mathbf{R}^N)$| and |$\varepsilon \in (0,1]$|⁠. Proof of Lemma2.5. See Appendix B. Since |$\nu ^r$| in (2.26) satisfies |$\| \nu ^r \|=j+k$| we can apply Lemma 2.5 to the expectation in (2.26). Then the expectation involving single iterated stochastic integrals of the length no less than |$3$| appearing in |${\prod _{l=1}^k} Y_1^{\beta _l,\alpha _l}$| of Lemma 2.4 will be at least |$O(\varepsilon ^6)$| if we use the regularity of |$\varphi $|⁠. In particular, the length of iterated stochastic integrals appearing in |$\prod _{l=1}^k Y_1^{\beta _l,\alpha _l}$| is greater than or equal to |$3$| whenever |$\sum _{l=1}^k \beta _l =3, k=2,3 $| or |$\sum _{l=1}^k \beta _l \geqslant 4 $|⁠. For |$\sum _{l=1}^k \beta _l = j$| with |$j=2,3, \ k=1$| the term |$\prod _{l=1}^k Y_1^{\beta _l,\alpha _l}$| includes double stochastic integrals such as |$I_{(j,0)}(1)$|⁠, |$I_{(0,j)}(1)$| (⁠|$\,j=1,\cdots ,d$|⁠), |$I_{(0,0)}(1)$| and iterated integrals of the length greater than |$2$|⁠, i.e. |$I_{J}(1)$|⁠, |$3\leqslant |J|\leqslant 4$|⁠, |$3 \leqslant \| J \| \leqslant 4$|⁠, |$J \in \{0,1,\cdots ,d \}^{|J|}$|⁠. Then, for the case |$j=2,3$|⁠, |$k=1$|⁠, the length of iterated stochastic integrals in |$\prod _{l=1}^k Y_1^{\beta _l,\alpha _l}$| will be greater than or equal to |$3$| except for |$I_{(j,0)}(1)$|⁠, |$I_{(0,j)}(1)$|⁠, |$ I_{(0,0)}(1)$|⁠. In addition, when |$\sum _{l=1}^k \beta _l =2,\ k=2$|⁠, the length of iterated stochastic integrals appearing in |$\prod _{l=1}^k Y_1^{\beta _l,\alpha _l}$| is no less than |$3$| except for |$ I_{(0,0)}(1)$| since we have $$\begin{eqnarray*} && I_{(j_1, j_2)}(1) I_{(j_3, j_4)}(1) \\ & = & I_{(j_1, j_2,j_3,j_4)}(1) + I_{(j_3, j_1, j_2, j_4)}(1) + I_{(j_1, j_3, j_2, j_4)}(1) + I_{(j_3, j_4, j_1, j_2)}(1) + I_{(j_1, j_3, j_4, j_2)}(1) + I_{(j_3, j_1, j_4, j_2)}(1) \nonumber \\ && + I_{(0, j_2, j_4)}(1) {\textbf{1}}_{j_1 =j_3 \neq 0} + I_{(j_1, 0, j_4)}(1) {\textbf{1}}_{j_2 =j_3 \neq 0} + I_{(0, j_4, j_2)}(1) {\textbf{1}}_{j_1 =j_3 \neq 0} \nonumber \\ && + I_{(j_3, 0, j_2)}(1) {\textbf{1}}_{j_1 =j_4 \neq 0} + I_{(j_1, j_3, 0)}(1) {\textbf{1}}_{j_2 =j_4 \neq 0} + I_{(j_3, j_1, 0)}(1) {\textbf{1}}_{j_2 =j_4 \neq 0} \nonumber \\ && + I_{(0,0)}(1) {\textbf{1}}_{j_1 = j_3 \neq 0, j_2 = j_4 \neq 0}, \nonumber \end{eqnarray*}$$ for |$j_1,\cdots ,j_4 \in \{0,1,\cdots ,d \}$|⁠. In summary, we have the following simple |$O(\varepsilon ^6)$|-small time expansion. Lemma 2.6 It holds that $$\begin{align} &\; E [\varphi ({X}_1^{x,\varepsilon})] = E [\varphi(\bar{X}_1^{x,\varepsilon})] + \sum_{i=1}^K E \left[ \partial_{i} \varphi ( \bar{X}_1^{x,\varepsilon} ) \sum_{0 \leqslant j_1,j_2 \leqslant d} \varepsilon^{\| (j_1,j_2)\|} \hat{V}_{j_1} V_{j_2}^{i}(x) I_{(j_1,j_2)}(1) \right] \end{align}$$(2.32) $$\begin{align} &\qquad +\frac{1}{2} \varepsilon^4 \sum_{i_1,i_2 =1}^K E \left[ \partial_{i_1} \partial_{i_2} \varphi ( \bar{X}_1^{x,\varepsilon} ) \sum_{1 \leqslant j_1, j_2 \leqslant d} \hat{V}_{j_1} V_{j_2}^{i_1}(x) \hat{V}_{j_1} V_{j_2}^{i_2}(x)I_{(0,0)}(1)\right] \\ &\qquad + \mathcal{E}_1^{\varphi,\varepsilon} (x)+ \mathcal{E}_2^{\varphi,\varepsilon} (x), \nonumber \end{align}$$(2.33) where |$\mathcal{E}_1^{\varphi ,\varepsilon } (x)$| satisfies $$\begin{equation*} \sup_{x \in \mathbf{R}^N} | \mathcal{E}_1^{\varphi,\varepsilon}(x) | \leqslant C \| \varphi \|_{\infty}\varepsilon^6\end{equation*}$$ for some |$C>0$|⁠, which are independent of |$\varphi $| and |$t$| by Lemma 2.4 and |$\mathcal{E}_2^{\varphi ,\varepsilon } (x)$| satisfies $$\begin{equation*}\mathcal{E}_2^{\varphi,\varepsilon}(x)=\varepsilon^6 \sum_{\lambda_l, l \leqslant m} E[ \partial^{\lambda_l} \varphi (\bar{X}_1^{x,\varepsilon})] g_{\lambda_l}(\varepsilon,x)\end{equation*}$$ for some |$m\geqslant 1$|⁠, |$\lambda _l \in \{1,\cdots ,N \}^{\,j}$|⁠, |$j\geqslant 1$|⁠, and the function |$g_{\lambda _l}$| is such that |$g_{\lambda _l}(\cdot ,x)$| is nondecreasing for any fixed |$x \in \mathbf{R}^N$| and |$g_{\lambda _l}(\varepsilon ,\cdot )$| is bounded for any fixed |$\varepsilon \in (0,1]$|⁠. Using integration by parts and a technique of distributions on Wiener space with |$X_1^{x,\sqrt{t}} \ {\overset{\mathrm{Law}}{\sim }} X_t^{x}$|⁠, we obtain simple representations for the terms in (2.32) and (2.33) using Malliavin weights, which are given as a polynomial of Gaussian random variables. Lemma 2.7 We have $$\begin{eqnarray*} &&\sum_{i=1}^K E \Bigg[ \partial_{i} \varphi ( \bar{X}_1^{x,\sqrt{t}} ) \sum_{0 \leqslant j_1,j_2 \leqslant d} {t}^{\frac{1}{2} \| (j_1,j_2)\|} \hat{V}_{j_1} V_{j_2}^{i}(x) I_{(j_1,j_2)}(1) \Bigg]\\ &=&E \Bigg[ \varphi (\bar{X}_t^x ) \sum_{i=1}^K \sum_{0 \leqslant j_1,j_2 \leqslant d} \hat{V}_{j_1} V_{j_2}^{i}(x) H_{(i)} \Bigg( \bar{X}_t ^x, \frac{1}{2}W_t^{j_1} W_t^{j_2} - \frac{t}{2} \mathbf{1}_{ j_1 = j_2 \neq 0} \Bigg) \nonumber \\ && \ \ \ + \sum_{i=1}^K \sum_{\alpha \in I} \sum_{\substack{ 0 \leqslant j \leqslant d \\ j \neq \alpha}} [V_{\alpha}, V_{j}]^i(x) H_{(i)} \Bigg( \bar{X}_t^x, \frac{1}{t} I_{(\alpha,0)}(t) W_t^{j} - \frac{1}{2} W_t^{\alpha} W_t^{j} \Bigg)\Bigg], \nonumber \end{eqnarray*}$$ and $$\begin{eqnarray*} &&\frac{1}{2} t^2 \sum_{i_1,i_2 =1}^K E \Bigg[ \partial_{i_1} \partial_{i_2} \varphi ( \bar{X}_1^{x,\sqrt{t}} ) \sum_{1 \leqslant j_1, j_2 \leqslant d} \hat{V}_{j_1} V_{j_2}^{i_1}(x) \hat{V}_{j_1} V_{j_2}^{i_2}(x)I_{(0,0)}(1)\Bigg]\\ &=& E \Bigg[ \varphi (\bar{X}_t^x ) \sum_{i_1,i_2 =1}^K \sum_{1 \leqslant j_1,j_2 \leqslant d} \hat{V}_{j_1} V_{j_2}^{i_1}(x) \hat{V}_{j_1} V_{j_2}^{i_2}(x) H_{(i_1,i_2)}( \bar{X}_t^x, 1) \frac{t^2}{4} \Bigg]. \end{eqnarray*}$$ Proof of Lemma2.7. See Appendix C. Finally, by applying Lemma 2.7 to (2.32) and (2.33) in Lemma 2.3 with the scaling |$\varepsilon =\sqrt{t}$|⁠, we easily obtain the result of Theorem 2.3. 2.4 Global approximation (main result) Hereafter, we assume |$T\geqslant 1$| and |$f: \mathbf{R}^N \rightarrow \mathbf{R}$| is bounded and measurable. Now we connect the local approximation operator |$n$|-times as |$(Q_{T/n})^n f (x),\; x \in \mathbf{R}^N,$| where |$(Q_t f)(\cdot )=E[\,f(\bar{X}_t^\cdot )\pi _t^\cdot ]$|⁠, in order to get the global approximation for |$P_Tf$|⁠. The following is our main result of the paper. Theorem 2.8 (Weak approximation). There exists a |$C>0$| such that $$\begin{eqnarray*} \Big\| P_T f - (Q_{T/n})^n f \Big\|_{\infty} \leqslant C \|\,f \|_{\infty} \frac{1}{n^2}, \end{eqnarray*}$$ for any bounded measurable function |$f:\mathbf{R}^{N} \rightarrow \mathbf{R}$| and |$n\geqslant 1$|⁠. 2.4.1 Proof of Theorem 2.8 Since the difference |$P_T f(x) - (Q_{T/n})^n f (x)$| is decomposed as $$\begin{align*} P_T f(x) - (Q_{T/n})^n f (x) = \sum_{k=0}^{n-1} (Q_{T/n})^{k} (P_{T/n} - Q_{T/n}) P_{T- {(k+1)T/n}} f(x), \ \ \ x \in \mathbf{R}^N, \end{align*}$$ it suffices to show $$\begin{eqnarray} \| (Q_{T/n})^{k} (P_{T/n} - Q_{T/n}) P_{T- {(k+1)T/n}} f \|_{\infty}=O \Big(\frac{1}{n^{3}}\Big) \end{eqnarray}$$(2.34) to prove the scheme |$(Q_{T/n})^n f(x)$| is the second order, i.e. $$\begin{eqnarray} \| P_T f - (Q_{T/n})^n f \|_{\infty} = O \Big(\frac{1}{n^{3}} \Big) \times n = O \Big(\frac{1}{n^{2}}\Big). \end{eqnarray}$$(2.35) Note that when |$f$| is sufficiently smooth (e.g. |$f \in C_b^\infty (\mathbf{R}^N)$|⁠) the approximation (2.35) is easily obtained using the following estimate for (2.34): $$\begin{eqnarray} \| (Q_{T/n})^{k} (P_{T/n} - Q_{T/n}) P_{T- {(k+1)T/n}} f\|_{\infty} \leqslant C \sum_{i=0}^m \| \nabla^i f \|_{\infty} \frac{1}{n^3} =O \Big( \frac{1}{n^3} \Big), \end{eqnarray}$$(2.36) which follows from the estimate |$\| Q_{T/n} g \|_{\infty } \leqslant C \| g \|_{\infty }$| for a bounded function |$g$| and from the local error decomposition $$\begin{eqnarray*} (P_{T/n} - Q_{T/n}) P_{T- {(k+1)T/n}} f(x)= \mathcal{E}_1^{P_{T- {(k+1)T/n}} f} (T/n,x)+ \mathcal{E}_2^{P_{T- {(k+1)T/n}} f} (T/n,x), \ \ \ x\in \mathbb{R}^N, \end{eqnarray*}$$ in Theorem 2.3 with the gradient estimate of |$P_{T- {(k+1)T/n}} f$| with the derivatives of |$f$|⁠: $$\begin{eqnarray} \| \partial^\beta P_{T- {(k+1)T/n}} f \|_{\infty} \leqslant C \sum_{i=1}^\ell \| \nabla^i f \|_{\infty}. \end{eqnarray}$$(2.37) Since the estimate (2.36) strongly depends on the smoothness of |$f$|⁠, it cannot be applied directly for the case that |$f$| is only bounded and measurable (i.e. nonsmooth). In fact, when |$f$| is bounded measurable, (2.37) does not hold and the gradient bound will be $$\begin{eqnarray} \| \partial^{\beta} P_{T- {(k+1)T/n}} f \|_{\infty} \leqslant C \|\,f \|_{\infty} \frac{1}{(T-(k+1)T/n)^{\nu}} \end{eqnarray}$$(2.38) for some |$\nu > 0$| by Kusuoka & Stroock (1985). Then the upper bound in (2.38) becomes large when |$k$| goes near |$n$|⁠, which suggests that some different approach is needed to show (2.34) and (2.35) under the nonsmooth condition on |$f$|⁠, unlike the smooth test function case. So, from now on, we take another strategy to show (2.34) and (2.35) under the condition that |$f$| is bounded and measurable. Note that |$(Q_{T/n})^n f(x)$| is given by $$\begin{align} (Q_{T/n})^n f(x)=E\left[\,f( \bar{X}_T^{(n),x} ) \prod_{k=1}^n {\pi}_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta_{kT/n}-\theta_{(k-1)T/n}) \right], \end{align}$$(2.39) where |${\bar{X}_{jT/n}^{(n),x}=\left (\bar{X}_{jT/n}^{(n),x,1},\ldots ,\bar{X}_{jT/n}^{(n),x,N}\right )}$|⁠, |$j=1,\cdots ,n$| is a Markov chain defined by $$\begin{align} \bar{X}_{jT/n}^{(n),x,i_1} =&\; \bar{X}_{(j-1)T/n}^{(n),x, i_1} + V_0^{i_1} \left( \bar{X}_{(j-1)T/n}^{(n),x} \right)\frac{T}{n} + \sum_{l=1}^d V_l^{i_1} \left( \bar{X}_{(j-1)T/n}^{(n),x} \right) \left( W_{jT/n}^l - W_{(j-1)T/n}^l\right),\; i_1 =1, \ldots, K, \end{align}$$(2.40) $$\begin{align} \bar{X}_{jT/n}^{(n),x, i_2}=&\; \bar{X}_{(j-1)T/n}^{(n),x,i_2} + V_0^{i_2} \left( \bar{X}_{(j-1)T/n}^{(n),x} \right)\frac{T}{n} +\sum_{l=0}^d \hat{V}_l V_0^{i_2} \left( \bar{X}_{(j-1)T/n}^{(n),x} \right) \int_{(j-1)T/n}^{jT/n} W_s^l \textrm{d}s, \; i_2 = K+1, \ldots, N, \end{align}$$(2.41) and |${\pi }_t^x(\theta _t)$| denotes |${\pi }_t^x$| with the dependence of $$\begin{align} \theta_t = \left( W_t^1, \ldots, W_t^d, \int_0^t W_s^{\alpha_1} \textrm{d}s, \ldots, \ \int_0^t W_s^{\alpha_l} \textrm{d}s \right), \end{align}$$(2.42) where |$\alpha _1, \ldots , \alpha _l$| are all elements in |$I:= {\bigcup _{m=K+1}^N} \{ \alpha \in \{ 1, 2,\ldots , d \}: \ [V_{\alpha }, V_0 ]^m (\cdot ) \neq 0 \}$|⁠. Then we have $$\begin{align*} &(Q_{T/n})^{k} (P_{T/n} - Q_{T/n}) P_{T- {(k+1)T/n}} f(x) \\ = &\; (Q_{T/n})^{k} \{ \mathcal{E}_1^{P_{T- {(k+1)T/n}} f}(T/n,x) + \mathcal{E}_2^{P_{T- {(k+1)T/n}} f} ( T/n,x) \} \nonumber \\ =&E\left[ \mathcal{E}_1^{P_{T- {(k+1)T/n}} f}(T/n,\bar{X}_{kT/n}^{(n),x}) \prod_{j=1}^k \pi_{T/n}^{ \bar{X}_{(j-1)T/n}^{(n),x}} \right] +E\left[ \mathcal{E}_2^{P_{T- {(k+1)T/n}} f}(T/n,\bar{X}_{kT/n}^{(n),x}) \prod_{j=1}^k \pi_{T/n}^{ \bar{X}_{(j-1)T/n}^{(n),x}} \right]. \end{align*}$$ By Lemma 2.4 we easily see $$\begin{align} \sup_{x \in \mathbf{R}^N} \left| E\left[ \mathcal{E}_1^{P_{T- {(k+1)T/n}} f}(T/n,\bar{X}_{kT/n}^{(n),x}) \prod_{j=1}^k \pi_{T/n}^{ \bar{X}_{(j-1)T/n}^{(n),x}}\right]\right| &\leqslant C \| P_{T- {(k+1)T/n}} f \|_\infty \frac{1}{n^3} \leqslant C \|\,f \|_\infty \frac{1}{n^3}. \end{align}$$(2.43) Also we have $$\begin{align*} &E\left[ \mathcal{E}_2^{P_{T- {(k+1)T/n}} f}(T/n,\bar{X}_{kT/n}^{(n),x}) \prod_{j=1}^k \pi_{T/n}^{ \bar{X}_{(j-1)T/n}^{(n),x}} \right]\\ =& \left(\frac{T}{n}\right)^3 \sum_{\lambda_l,l\leqslant m} E\left[ \partial^{\lambda_l} P_{T- {(k+1)T/n}} f(\bar{X}_{(k+1)T/n}^{(n),x}) {\cal G}_{\lambda_l}(T/n,\bar{X}_{kT/n}^{(n),x}) \prod_{j=1}^k \pi_{T/n}^{ \bar{X}_{(j-1)T/n}^{(n),x}}\right]. \end{align*}$$ Then it suffices to show that for any multi-index |$\beta $| and |$G \in \mathbf{D}^{\infty }$| such that for all |$l\geqslant 1$|⁠, |$p<\infty $|⁠, |$\| G \|_{l,p} \leqslant c$| for some |$c>0$|⁠, there is a |$C>0$| such that $$\begin{eqnarray} \sup_{x \in \mathbf{R}^N} | E[ \partial^{\beta} P_{T- {(k+1)T/n}}\ f(\bar{X}_{(k+1)T/n}^{(n),x}) G] | \leqslant C \|\,f \|_{\infty}, \end{eqnarray}$$(2.44) for all |$n\geqslant 1$| and |$k=0,1,\cdots ,n-1$|⁠. Recall again that for a multi-index |$\beta $| we can choose |$C>0$| and |$\nu>0$|⁠, $$\begin{eqnarray} \| \partial^{\beta} P_{T- s}\ f \|_{\infty} \leqslant C \|\,f \|_{\infty} \frac{1}{(T-s)^{\nu}}, \end{eqnarray}$$(2.45) for all |$s \in (0,T)$|⁠, by Kusuoka & Stroock (1985) under the hypoelliptic condition. Then the estimate (2.44) is immediately obtained when |$k=0,1,\cdots ,[n/2]-1$| (i.e. |$T/2 < T-s \leqslant T$|⁠). However, when |$k=[n/2],\cdots ,n-1$|⁠, (2.45) cannot be applied to the left-hand side of (2.44) because the upper bound of (2.45) may be huge according to the values of |$k$| and |$n$|⁠. Although we have this problem we can still show (2.44) by a truncation argument. In order to prove (2.44) we first prepare the following lemma. Lemma 2.9 For |$k=[n/2]+1,\cdots ,n$|⁠, |$\bar{X}_{kT/n}^{(n),x}$| and |${X}_{kT/n}^{x}$| are nondegenerate in Malliavin's sense. In particular, for any |$p>1$| there exist a constant |$C_1>0$| and an integer |$m$| such that $$\begin{align} \sup_{k=[n/2]+1,\cdots,n} E [ ( \det \sigma^{\bar{X}_{kT/n}^{(n),x}})^{-p}] \leqslant C_1 n^{m}, \end{align}$$(2.46) for all |$x \in \mathbf{R}^N$|⁠. Furthermore, for any |$p>1$| there exist |$C_2>0$| and |$\nu>0$| such that $$\begin{align} \sup_{T/2\leqslant t \leqslant T}E [ ( \det \sigma^{{X}_{t}^x})^{-p}] \leqslant C_2 \frac{1}{T^\nu}, \end{align}$$(2.47) for all |$x \in \mathbf{R}^N$|⁠. Proof of Lemma2.9. See Appendix D. While the functional |${\bar{X}_{kT/n}^{(n),x}}$|⁠, |$k=[n/2]+1,\ldots ,n$|⁠, does not degenerate at least the moment of the inverse of the Malliavin covariance depends on the number of discretization |$n$| and may be huge. However, we can use the bound (2.47) in order to get (2.44) on the set that |${\bar{X}_{kT/n}^{(n),x}}$| is close to the original process |${X}_{kT/n}^x$|⁠. Furthermore, the probability of the complement of the set |${\bar{X}_{kT/n}^{(n),x}}$| is close to |${X}_{kT/n}^x$| will be very small and negligible. So we introduce a truncation functional $$\begin{equation*} r_{kT/n}:= {\det \sigma^{\bar{X}_{kT/n}^{(n),x}}/\det \sigma^{X_{kT/n}^x} \ -1}.\end{equation*}$$ Let |$\phi \in C_b^\infty (\mathbf{R})$| be a function such that |$\phi (x)=1$| for |$|x| \leqslant 1/4$|⁠, |$\phi (x)=0$| for |$|x|\geqslant 1/2$| and |$\phi (x) \in (0,1)$| for |$|x|\in (1/4,1/2)$|⁠. Let |$\psi _{kT/n}$| be a truncation given by |$\psi _{kT/n}=\phi (r_{kT/n}) \in \mathbf{D}^{\infty }$|⁠. We have the following result, which leads to the proof of (2.44). Lemma 2.10 Let |$\beta \in \{1, \ldots , N\}^e, e \in \mathbf{N}$|⁠, be a multi-index, and let |$\nu>1$| and |$G \in \mathbf{D}^{\infty }$| such that for all |$l \geqslant 1$|⁠, |$p<\infty $|⁠, |$\| G \|_{l,p} \leqslant c$| for some |$c>0$|⁠. Then, there exist |$C_1,C_2>0$| such that $$\begin{eqnarray} &&\sup_{x \in \mathbf{R}^N}| E[ \partial^{\beta} P_{T-(k+1)T/n} f( \bar{X}_{(k+1)T/n}^{(n),x} ) G \psi_{(k+1)T/n} ] | \leqslant C_1 \|\,f \|_{\infty}, \end{eqnarray}$$(2.48) $$\begin{eqnarray} &&\sup_{x \in \mathbf{R}^N}| E[ \partial^{\beta} P_{T-(k+1)T/n} f( \bar{X}_{(k+1)T/n}^{(n),x} ) G (1-\psi_{(k+1)T/n}) ] | \leqslant C_2 \frac{1}{n^\nu} \|\,f \|_{\infty}, \end{eqnarray}$$(2.49) for all |$n\geqslant 1$|⁠, |$k=[n/2],\cdots ,n-1$|⁠, and a bounded measurable function |$f: \mathbf{R}^N \to \mathbf{R}$|⁠. Proof of Lemma 2.10. By Lemma 2.9 and (2.4), we have $$\begin{align*} {}_{-\infty} \langle \partial^{\alpha} S(\bar{X}_{(k+1)T/n}^{(n),x}), \eta \rangle_{\infty} = {}_{-\infty} \langle S(\bar{X}_{(k+1)T/n}^{(n),x}), H_\alpha (\bar{X}_{(k+1)T/n}^{(n),x},\eta) \rangle_{\infty} \end{align*}$$ for all |$k=[n/2], \ldots , n-1$|⁠, multi-index |$\alpha $|⁠, |$S \in \mathcal{S}^{\prime}(\mathbf{R}^N)$| and |$\eta \in \mathbf{D}^{\infty }$|⁠. Then it holds that $$\begin{align*} &| E[ \partial^{\beta} P_{T-(k+1)T/n} f( \bar{X}_{(k+1)T/n}^{(n),x} ) G \psi_{(k+1)T/n} ] | \\ = &\; \left| E \Big[ P_{T-(k+1)T/n} f( \bar{X}_{(k+1)T/n}^{(n),x} ) H_{\beta} ( \bar{X}_{(k+1)T/n}^{(n),x}, G \psi_{(k+1)T/n} ) \mathbf{1}_{ \{ | r_{(k+1)T/n} | \leqslant 1/2 \}} \Big] \right| \\ \leqslant &\; K_1 \|\,f \|_{\infty} E \left[ \Big| H_{\beta} ( \bar{X}_{(k+1)T/n}^{(n),x}, G \psi_{(k+1)T/n} ) \mathbf{1}_{ \{ | r_{(k+1)T/n} | \leqslant 1/2 \}} \Big| \right] \\ \leqslant &\; K_2 \|\,f \|_{\infty} \Big\| ( \det \sigma^{\bar{X}_{(k+1)T/n}^{(n),x}})^{-1} {\textbf{1}}_{\big\{ {\det \sigma^{ \bar{X}_{(k+1)T/n}^{(n),x}} \geqslant{1}/{2} \det \sigma^{{X}_{(k+1)T/n}^{x}} \big\}}} \Big\|_p^m\\ \leqslant &\; K_3 \|\,f \|_{\infty} \| ( \det \sigma^{X_{(k+1)T/n}^x})^{-1} \|_p^m, \end{align*}$$ for some |$K_1,K_2,K_3>0$|⁠, |$p>1$| and an integer |$m$|⁠, which are independent of |$x$|⁠. Then by Kusuoka–Stroock’s bound (2.47) we can choose a constant |$C_1>0$| in (2.48), which is independent of |$n$|⁠. Furthermore, again by Lemma 2.9 and especially by (2.46), we have $$\begin{align*} &| E[ \partial^{\beta} P_{T-(k+1)T/n} f ( \bar{X}_{(k+1)T/n}^{(n),x} ) G (1-\psi_{(k+1)T/n}) ] | \\ =&\; | E[ \partial^{\beta} P_{T-(k+1)T/n}f ( \bar{X}_{(k+1)T/n}^{(n),x} ) G (1-\psi_{(k+1)T/n}) \mathbf{1}_{ \{ | r_{(k+1)T/n} | \geqslant 1/4 \}} ] | \\ \leqslant&\; K_4 \|\,f \|_{ \infty} n^{m} \mathbb{P} \left( \Big| (\det \sigma^{{X}^{x}_{(k+1)T/n}} )^{-1} \det \sigma^{\bar{X}^{(n),x}_{(k+1)T/n}} -1 \Big| \geqslant 1/4 \right), \end{align*}$$ for a constant |$K_4>0$| and an integer |$m$|⁠, which are independent of |$x$|⁠. Note that it holds $$\begin{eqnarray} &&\mathbb{P} \left( \Big| (\det \sigma^{{X}^{x}_{(k+1)T/n}} )^{-1} \det \sigma^{\bar{X}^{(n),x}_{(k+1)T/n}} -1 \Big| \geqslant 1/4 \right) \nonumber \\ &=& \mathbb{P} \left( \Big| \det \sigma^{\bar{X}^{(n),x}_{(k+1)T/n}} -\det \sigma^{{X}_{(k+1)T/n}^{x}} \Big| \geqslant (\det \sigma^{X^x_{(k+1)T/n}} )/4 \right)\nonumber\\ &\leqslant& \mathbb{P} \left( ( \det \sigma^{{X}^x_{(k+1)T/n}} )^{-1} \geqslant n^{1/4} \right)+\mathbb{P} \left( \Big| \det \sigma^{\bar{X}^{(n),x}_{(k+1)T/n}} -\det \sigma^{{X}^x_{(k+1)T/n}} \Big| \geqslant 1/(4n^{1/4}) \right)\nonumber\\ &\leqslant& n^{-r/4} E [ |\! \det \sigma^{{X}^x_{(k+1)T/n}} |^{-r} ] + (4n^{1/4})^r E [ |\! \det \sigma^{\bar{X}^{(n),x}_{(k+1)T/n}} -\det \sigma^{{X}^x_{(k+1)T/n}} |^{r} ] \end{eqnarray}$$(2.50) for all |$r>1$|⁠. Also, using the similar argument of Bally & Talay (1996), we get for all |$r>1$|⁠, $$\begin{eqnarray} E [ |\!\det \sigma^{\bar{X}^{(n),x}_{(k+1)T/n}} -\det \sigma^{{X}^x_{(k+1)T/n}} |^{r} ] = O( n^{-r/2} ). \end{eqnarray}$$(2.51) Then by applying Kusuoka–Stroock’s bound (2.47) and (2.51) to (2.50), we obtain $$\begin{eqnarray*} \mathbb{P} \left( \Big| (\det \sigma^{{X}_{(k+1)T/n}^x} )^{-1} \det \sigma^{\bar{X}^{(n),x}_{(k+1)T/n}} -1 \Big| \geqslant 1/4 \right) = O( n^{-\eta}), \ \ \textrm{for all}\ \eta>1. \end{eqnarray*}$$ Therefore, we have $$\begin{eqnarray*} &&\sup_{x\in \mathbf{R}^N} \Big| {E} [\partial^{\beta} P_{T-(k+1)T/n}f ( \bar{X}^{(n),x}_{(k+1)T/n} ) G ( 1 - \psi_{(k+1)T/n})] \Big| \leqslant K_5 \|\,f \|_{\infty} n^{ m - \eta}, \textrm{for all}\ \eta>1, \end{eqnarray*}$$ for a constant |$K_5>0$|⁠. Since we can choose |$\eta $| arbitrarily we obtain (2.49). |$\Box $| Now, we return to the proof of Theorem 2.8. Applying Lemma 2.10 to $$\begin{eqnarray*} &&E[ \partial^{\beta} P_{T- {(k+1)T/n}} f(\bar{X}_{(k+1)T/n}^{(n),x}) G]\big|\\ &=&E[ \partial^{\beta} P_{T- {(k+1)T/n}} f(\bar{X}_{(k+1)T/n}^{(n),x}) G(1-\psi_{(k+1)T/n})]+E[ \partial^{\beta} P_{T- {(k+1)T/n}} f(\bar{X}_{(k+1)T/n}^{(n),x}) G \psi_{(k+1)T/n}], \end{eqnarray*}$$ with |$G={\cal G}_{\beta }(T/n,\bar{X}_{kT/n}^{(n),x}) \prod _{j=1}^k \pi _{T/n}^{ \bar{X}_{(j-1)T/n}^{(n),x}}$|⁠, we easily obtain (2.44) and $$\begin{align*} &\Bigg| E\Bigg[ \mathcal{E}_2^{P_{T- {(k+1)T/n}} f}(T/n,\bar{X}_{kT/n}^{(n),x}) \prod_{j=1}^k \pi_{T/n}^{ \bar{X}_{(j-1)T/n}^{(n),x}} \Bigg]\Bigg|\\ \leqslant & \left(\frac{T}{n}\right)^3 \sum_{\lambda_l,l\leqslant m} \Bigg| E\Bigg[ \partial^{\lambda_l} P_{T- {(k+1)T/n}} f(\bar{X}_{(k+1)T/n}^{(n),x}) {\cal G}_{\lambda_l}(T/n,\bar{X}_{kT/n}^{(n),x}) \prod_{j=1}^k \pi_{T/n}^{ \bar{X}_{(j-1)T/n}^{(n),x}}\Bigg]\Bigg| \\ \leqslant & C \|\,f \|_\infty \frac{1}{n^3}, \end{align*}$$ for all |$x \in \mathbf{R}^N$|⁠. Since we have already seen (in (2.43)): $$\begin{eqnarray*} \sup_{x \in \mathbf{R}^N} \Bigg| E\Bigg[ \mathcal{E}_1^{P_{T- {(k+1)T/n}} f}(T/n,\bar{X}_{kT/n}^{(n),x}) \prod_{j=1}^k \pi_{T/n}^{ \bar{X}_{(j-1)T/n}^{(n),x}}\Bigg]\Bigg| \leqslant C \|\,f \|_\infty \frac{1}{n^3}, \end{eqnarray*}$$ we obtain $$\begin{align*} &|(Q_{T/n})^{k} (P_{T/n} - Q_{T/n}) P_{T- {(k+1)T/n}} f(x) | \\ \leqslant &\Bigg|E\Bigg[ \mathcal{E}_1^{P_{T- {(k+1)T/n}} f}(T/n,\bar{X}_{kT/n}^{(n),x}) \prod_{j=1}^k \pi_{T/n}^{ \bar{X}_{(j-1)T/n}^{(n),x}} \Bigg] \Bigg| +\Bigg|E\Bigg[ \mathcal{E}_2^{P_{T- {(k+1)T/n}} f}(T/n,\bar{X}_{kT/n}^{(n),x}) \prod_{j=1}^k \pi_{T/n}^{ \bar{X}_{(j-1)T/n}^{(n),x}} \Bigg]\Bigg| \\ \leqslant & C \|\,f \|_\infty \frac{1}{n^3}, \end{align*}$$ for all |$x \in \mathbf{R}^N$|⁠. Therefore, we finally have $$\begin{eqnarray*} \| P_T f(x) - (Q_{T/n})^n f \|_{\infty} &\leqslant& \sum_{k=0}^{n-1} \| (Q_{T/n})^{k} (P_{T/n} - Q_{T/n}) P_{T- {(k+1)T/n}} f \|_{\infty} \\ & \leqslant & C\|\,f \|_{\infty} \sum_{k=0}^{n-1} \frac{1}{n^3}= C \|\,f \|_{\infty} \frac{1}{n^2}. \end{eqnarray*}$$ 2.4.2 Some extension The bounded measurable condition on the test function is replaced with the Lipschitz continuous one. In such a case the statement of Theorem 2.8 will be as follows. Corollary 2.11 There exists a |$C>0$| such that $$\begin{eqnarray*} \Big\| P_T f- (Q_{T/n})^n f \Big\|_{\infty} \leqslant C C_{\mathrm{Lip}}[\,f] \frac{1}{n^2}, \end{eqnarray*}$$ for any Lipschitz continuous function |$f:\mathbf{R}^{N} \rightarrow \mathbf{R}$| and |$n\geqslant 1$|⁠, where |$C_{\mathrm{Lip}}[\,f]$| denotes the Lipschitz constant of |$f$|⁠. Proof of Corollary 2.11. Let |$\beta $| be a multi-index and |$G \in \mathbf{D}^{\infty }$| a Wiener functional such that for all |$l\geqslant 1$|⁠, |$p<\infty $|⁠, there exists a |$c>0$| such that |$\| G \|_{l,p} \leqslant c$|⁠. Then, there exists a |$C>0$| such that $$\begin{eqnarray} | E[ \partial^{\beta} P_{T- {(k+1)T/n}} f(\bar{X}_{(k+1)T/n}^{(n),x}) G] | \leqslant C C_{\mathrm{Lip}}[\,f], \end{eqnarray}$$(2.52) for any Lipschitz continuous function |$f:\mathbf{R}^{N} \rightarrow \mathbf{R}$|⁠, |$x \in \mathbf{R}^N$|⁠, |$n\geqslant 1$| and |$k=0,1,\cdots ,n-1$|⁠. The assertion follows by the same argument as in Theorem 2.8. 2.5 Exact formula of the hypoelliptic weight In this section we give the exact formula for the Malliavin weight |$\pi _t^x$| in the local approximation operator in Theorem 2.3 and comment on the scheme from the perspective of implementations. The weight |${\pi }_t^x$| is explicitly computed by a linear combination of polynomials of Gaussian random variables while it involves Skorohod integrals. This is because we are able to obtain explicit forms of the Malliavin derivative and the Malliavin covariance of |$\bar{X}_t^x$| appearing in the type of weight |$H_{(i)}( \bar{X}_t^x,G)$| with |$G$| given by a polynomial of |$W_t^k$| and |${\int _0^t W_s^l \textrm{d}s}$|⁠, |$1\leqslant k,l \leqslant d$|⁠. In fact we have the following formulas for the functionals in |${\pi }_t^x$| in Theorem 2.3. Proposition 2.12 (Exact formula for hypoelliptic weight). It holds that for |$1 \leqslant i, i_1, i_2 \leqslant K, \; 0 \leqslant j_1, j_2 \leqslant d $|⁠, $$\begin{align} & H_{(i)} \Big( \bar{X}_t^x, \frac{1}{2}{W_t^{j_1}W_t^{j_2}}- \frac{t}{2} \mathbf{1}_{ j_1 = j_2 \neq 0} \Big) \nonumber \\ = &\; \Big( \frac{1}{2}{W_t^{j_1}W_t^{j_2}}- \frac{t}{2} \mathbf{1}_{ j_1 = j_2 \neq 0} \Big) \sum_{m=1}^d \left\{ \sum_{l=1}^K \gamma^{\bar{X}_t^x}_{il} V_m^l (x) W_t^m + \sum_{l=K + 1}^N \gamma^{\bar{X}_t^x}_{il} \hat{V}_m V_0^l(x) \int_0^t W_s^m \textrm{d}s \right\} \nonumber\\ & - \frac{1}{2} \sum_{l=1}^K \gamma^{\bar{X}_t^x}_{il} \left( V_{j_2}^l(x) W_t^{j_1} + V_{j_1}^l(x) W_t^{j_2} \right) t -\frac{1}{4} \sum_{l=K + 1}^N \gamma^{\bar{X}_t^x}_{il} \left( \hat{V}_{j_2} V_0^l (x) W_t^{j_1} + \hat{V}_{j_1} V_0^l (x) W_t^{j_2} \right) t^2 \end{align}$$(2.53) and $$\begin{align} & H_{(i_1,i_2)} ( \bar{X}_t^x, 1) = H_{(i_1)} (\bar{X}_t^x,1 ) H_{(i_2)} (\bar{X}_t^x,1 ) - \gamma^{\bar{X}_t^x}_{i_1 i_2}, \end{align}$$(2.54) with $$\begin{align} H_{(i)} (\bar{X}_t^x,1 ) = \sum_{m=1}^d \left\{ \sum_{l=1}^K \gamma^{\bar{X}_t^x}_{il}V_m^l (x) W_t^m + \sum_{l=K+1}^N \gamma^{\bar{X}_t^x}_{il} \hat{V}_m V_0^l(x) \int_0^t W_s^m \textrm{d}s \right\}. \nonumber \end{align}$$ Also, for |$1 \leqslant i \leqslant K, \alpha \in I$| and |$0 \leqslant j \leqslant d, \; j \neq \alpha $|⁠, we have $$\begin{align} & H_{(i)} \left( \bar{X}_t^x, \frac{1}{t} \left(\int_0^t W_s^{\alpha} \textrm{d}s \right) W_t^{j}- \frac{1}{2} W_t^{\alpha} W_t^{\,j} \right) \nonumber\\ = &\; \Big( \frac{1}{t} \Big(\int_0^t W_s^{\alpha} \textrm{d}s \Big) W_t^{j}- \frac{1}{2} W_t^{\alpha} W_t^{\,j} \Big) \sum_{m=1}^d \left\{ \sum_{l=1}^K \gamma^{\bar{X}_t^x}_{i l} V_m^l (x) W_t^m + \sum_{l=K+1}^N \gamma^{\bar{X}_t^x}_{i l} \hat{V}_m V_0^l (x) \int_0^t W_s^m \textrm{d}s \right\} \nonumber\\ &\; - \left( \frac{1}{t} \Big(\int_0^t W_s^{\alpha} \textrm{d}s \Big) - \frac{1}{2} W_t^{\alpha} \right)\! \left\{ \sum_{l=1}^K \gamma^{\bar{X}_t^x}_{i l} V_j^l (x) t +\! \sum_{l= K+1}^N \gamma^{\bar{X}_t^x}_{i l} \frac{1}{2} \hat{V}_j V_0^l(x) t^2 \right\} -\! \sum_{l=K+1}^N \gamma^{\bar{X}_t^x}_{i l} \frac{1}{12} \hat{V}_\alpha V_0^l (x) W_t^{\,j} t^2. \end{align}$$(2.55) In particular, |$ \gamma ^{\bar{X}_t^x} = (\gamma ^{\bar{X}_t^x}_{i_1 i_2})$| is the deterministic |$N \times N$| symmetric matrix defined by |$\gamma ^{\bar{X}_t^x} = ( \sigma ^{\bar{X}_t^x} )^{-1}$|⁠, where |$\sigma ^{\bar{X}_t^x}$| is the |$N \times N$| symmetric positive definite matrix given by $$\begin{equation} \sigma^{\bar{X}_t^x}_{i_1 i_2} = \left\{ \begin{array}{lll} \sum_{l=1}^d V_l^{i_1} ( x ) V_l^{i_2} (x) t & 1 \leqslant i_1, i_2 \leqslant K \\ \sum_{l=1}^d V_l^{i_1} ( x ) \hat{V}_l V_0^{i_2} ( x )\frac{1}{2} t^2 & 1 \leqslant i_1 \leqslant K, \; K+1 \leqslant i_2 \leqslant N \\ \sum_{l=1}^d \hat{V}_l V_0^{i_1} ( x ) \hat{V}_l V_0^{i_2} (x ) \frac{1}{3} t^3 & K+1 \leqslant i_1, i_2 \leqslant N. \end{array} \right. \end{equation}$$(2.56) Proof of Proposition2.12. See Appendix E. The weight |${\pi }_t^x$| can be exactly simulated by either the Monte Carlo method or a quasi-Monte Carlo method. This is a key feature of our scheme. Also, it is different from the method of exact simulation of diffusions provided by Beskos & Roberts (2005) and the parametrix method proposed by Bally & Kohatsu-Higa (2015). Our scheme does not involve any implicit computation in the implementation. The algorithm is given as follows. Algorithm 2.13 (Monte Carlo method for the second-order scheme) Let |$\bar{X}_t^x(\theta _t)$| denote |$\bar{X}_t^x$| and recall that the stochastic process |$\{ \theta _t \}_{t \geqslant 0}$| is defined in (2.42). for|$j=1$|to|$M$|do |$ \bar{X}^{(n),x,[j]}_{0}=x$|⁠, |$\mathbb{W}^{[j]}=1$| for|$k=1$|to|$n$|do Simulate |$\theta ^{[j]}_{kT/n}-\theta ^{[j]}_{(k-1)T/n}$| Update |$\mathbb{W}^{[j]}=\mathbb{W}^{[j]} \times{\pi }_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x,[j]}}(\theta ^{[j]}_{kT/n}-\theta ^{[j]}_{(k-1)T/n})$| Compute |$ \bar{X}^{(n),x,[j]}_{kT/n}=\bar{X}_{T/n}^{\bar{X}^{(n),x,[j]}_{(k-1)T/n}}(\theta ^{[j]}_{kT/n}-\theta ^{[j]}_{(k-1)T/n})$| endfor endfor Return |$\frac{1}{M} \sum _{j=1}^M f(\bar{X}^{(n),x,[j]}_{T}) \mathbb{W}^{[j]}$|⁠. The Monte Carlo method provides an efficient scheme since the variance is finite as will be shown in Section 3.1. Also, some control variate methods are introduced in subsection Section 3.2 (3.2.1 and 3.2.2). The use of a quasi-Monte Carlo method is also recommended for some cases because it gives faster convergence. Numerical examples are given in Section 4 in order to check the effectiveness of the scheme. 3. Analysis of the numerical scheme In the following we show that the variance of the proposed second-order scheme is finite for any number of time steps. Moreover, in order to obtain faster convergence, we will introduce new control variate methods (in Section 3.2). 3.1 Variance estimation of the new scheme First, we mention the variance of the Euler–Maruyama scheme. The |$L^2$|-norm of the scheme is given as follows. Proposition 3.1 (The |$L^2$|-norm of the Euler–Maruyama scheme). Let |$\bar{X}_T^{\mathrm{EM},(n),x}$| be the |$n$|-time step Euler–Maruyama scheme of |${X}_T^{x}$| with the partition |$T/n$|⁠. There exist constants |$C,c>0$| such that $$\begin{align} \|\,f ( \bar{X}_T^{\mathrm{EM},(n),x} ) \|_{2} \leqslant C e^{c T}, \nonumber \end{align}$$ for any |$n \in \mathbf{N}$|⁠, |$x \in \mathbf{R}^N$| and measurable functions |$f: \mathbf{R}^N \rightarrow \mathbf{R}$|⁠, which have at most linear growth. Proof. See Proposition 7.2 of Pagès (2018). Next, we will show that the variance of the proposed scheme is finite for any discretization steps and has the same order to the above. Recall that the approximation |$(Q_{T/n})^n f(x)$| is defined as $$\begin{align} P_T f(x) \simeq (Q_{T/n})^n f(x)=E\left[\,f( \bar{X}_T^{(n),x} ) \prod_{k=1}^n {\pi}_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta_{kT/n}-\theta_{(k-1)T/n}) \right]. \end{align}$$(3.1) Then, the |$L^2$|-norm of the new scheme is given as follows. Proposition 3.2 (The |$L^2$|-norm of the proposed scheme). There exist constants |$C,c>0$| such that $$\begin{align} \left\|\,f( \bar{X}_T^{(n),x} ) \prod_{k=1}^n {\pi}_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta_{kT/n}-\theta_{(k-1)T/n}) \right\|_2 \leqslant C e^{c T}, \end{align}$$(3.2) for any |$n \in \mathbf{N}$|⁠, |$x \in \mathbf{R}^N$| and measurable functions |$f:\mathbf{R}^N \rightarrow \mathbf{R}$|⁠, which have at most linear growth. Proof. Using Hölder’s inequality, we obtain $$\begin{align} & \left\|\,f( \bar{X}_T^{(n),x} ) \prod_{k=1}^n {\pi}_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta_{kT/n}-\theta_{(k-1)T/n}) \right\|_2 \nonumber \\ \leqslant & \; \|\,f( \bar{X}_T^{(n),x} ) \|_{4} \left\| \prod_{k=1}^n {\pi}_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta_{kT/n}-\theta_{(k-1)T/n}) \right\|_{4}. \end{align}$$(3.3) We note that there exists a constant |$C>0$| such that |$E [ | \bar{X}_T^{(n),x} |^{4} ] \leqslant (1 +|x|^{4})e^{CT}$| due to a similar argument to give the estimation of the moments of the Euler–Maruyama scheme. Furthermore, there is a constant |$L>0$| such that |$|f(x)| \leqslant L (1 + |x|)$| for all |$x \in \mathbf{R}^N$| by the linear growth property of function |$f$|⁠. Hence, we have $$\begin{align} \|\,f ( \bar{X}_T^{(n),x} ) \|_{4} \leqslant L (1 + \| \bar{X}^{(n),x}_T \|_4 ) \leqslant C_1e^{C_2 T} \end{align}$$(3.4) for positive constants |$C_1$| and |$C_2$|⁠, which are uniform with respect to |$n$|⁠. Next, we will give the upper bound for $$\begin{align} \left\| \prod_{k=1}^n {\pi}_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta_{kT/n}-\theta_{(k-1)T/n}) \right\|_4. \end{align}$$(3.5) Let us denote the Brownian filtration by |$\{ \mathcal{F}_t\}_{t \in [0,T]}$|⁠. Then, by the tower property of conditional expectation, we have $$\begin{equation} E \left[ \left| \prod_{k=1}^n {\pi}_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta_{kT/n}-\theta_{(k-1)T/n}) \right|^4 \right] \end{equation}$$(3.6) $$\begin{align} = & E \left[ \left| \prod_{k=1}^{n-1} {\pi}_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta_{kT/n}-\theta_{(k-1)T/n}) \right|^4 \big. E \left[ \left| {\pi}_{T/n}^{\bar{X}_{(n-1)T/n}^{(n),x}}(\theta_{T}-\theta_{(n-1)T/n}) \right|^4 \bigg| \mathcal{F}_{(n-1)T/n} \right] \right]. \end{align}$$(3.7) Since |$\bar{X}_{(n-1)T/n}^{(n),x}$| is |$\mathcal{F}_{(n-1)T/n}$|-measurable and |$\theta _T - \theta _{(n-1)T/n}$| is independent of |$\mathcal{F}_{(n-1)T/n}$|⁠, we have $$\begin{align} \big. E \Big[ \big| {\pi}_{T/n}^{\bar{X}_{(n-1)T/n}^{(n),x}}(\theta_{T}-\theta_{(n-1)T/n}) \Big|^4 \big| \mathcal{F}_{(n-1)T/n} \Big] = \left. E \left[ \big| {\pi}_{T/n}^{z}(\theta_{T}-\theta_{(n-1)T/n}) \big|^4 \right] \right|{}_{z= \bar{X}_{(n-1)T/n}^{(n),x}}. \end{align}$$(3.8) Moreover, noting that $$\begin{align} E [ (W_t)^{n} ] = \begin{cases} 0 & ( n: \textrm{odd}) \\ \frac{ n!}{2^{n/2}(n/2)!}t^{n/2}& ( n: \textrm{even}), \end{cases} \end{align}$$(3.9) we are able to show that under our assumption there exists a constant |$C>0$| such that $$\begin{align} E \left[ \big| {\pi}_{T/n}^{z}(\theta_{T}-\theta_{(n-1)T/n}) \big|^4 \right] \leqslant 1 + C\frac{T}{n} \end{align}$$(3.10) for all |$z \in \mathbf{R}^N$| and hence, $$\begin{align} E \left[ \left| \prod_{k=1}^n {\pi}_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta_{kT/n}-\theta_{(k-1)T/n}) \right|^4 \right] \leqslant \left( 1 + C \frac{T}{n} \right) E \left[ \left| \prod_{k=1}^{n-1} {\pi}_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta_{kT/n}-\theta_{(k-1)T/n}) \right|^4 \right]. \nonumber \end{align}$$ Repeating the same procedure for the expectation appearing in the above right-hand side, we obtain $$\begin{align} E \left[ \left| \prod_{k=1}^n {\pi}_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta_{kT/n}-\theta_{(k-1)T/n}) \right|^4 \right] \leqslant \left( 1 + c \frac{T}{n} \right)^n \leqslant e^{c T} \end{align}$$(3.11) for |$c>0$|⁠. Then, we have $$\begin{align} \left\| \prod_{k=1}^n {\pi}_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta_{kT/n}-\theta_{(k-1)T/n}) \right\|_4 \leqslant e^{CT} \end{align}$$(3.12) for a positive constant |$C$|⁠. Applying the inequalities (3.4) and (3.12) to (3.3) we easily reach to the assertion. 3.2 Control variate methods for the new scheme Proposition 3.2 implies that an explosion of the variance will never occur in the proposed scheme as the number of time steps |$n \to \infty $|⁠, which is also checked in numerical examples in Section 4. While the variance is finite we can introduce some control variate methods in order to obtain a more efficient method. First, we quickly review the concept of control variate, which is one of the important variance reduction techniques. Now, we assume the situation that we aim at computing |$E[Y]$| by Monte Carlo simulation where |$Y$| is a random variable. Then, we choose another random variable |$Z$| and a constant |$b$| and define the control variate |$Y^b:= Y - b(Z - E[Z] )$|⁠, which is unbiased for the original random variable |$Y$|⁠. If it holds |$Var[Y] \geqslant Var[Y^b]$| it is preferable to compute |$E[Y^b]$| instead of |$E[Y]$|⁠. In fact it is well known that |$b^* = Cov[Y, Z] / Var[Z]$| minimizes |$Var[Y^b]$| and one has |$Var[Y^{b^*}] = (1-\rho ^2_{YZ}) Var[Y] \leqslant Var[Y]$|⁠, where |$\rho _{YZ}$| is the correlation between |$Y$| and |$Z$|⁠. Hence, in order to construct a desirable control variable, it is reasonable to choose a random variable |$Z$|⁠, which has a closed form of |$E[Z]$| and highly correlated with |$Y$|⁠. Based on the above discussion we introduce two control variate methods for the proposed numerical scheme (3.1). 3.2.1 Control variate 1 As the first control variate we consider $$\begin{align} Y_{(n)}^{CV1}:= f( \bar{X}_T^{(n),x} ) \prod_{k=1}^n {\pi}_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta_{kT/n}-\theta_{(k-1)T/n}) - \lambda \bigg( \prod_{k=1}^n {\pi}_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta_{kT/n}-\theta_{(k-1)T/n}) - 1 \bigg), \end{align}$$(3.13) where |$\lambda \in \mathbf{R}$| is defined by $$\begin{align} \lambda = \frac{Cov \Big[\,f( \bar{X}_T^{(n),x} ) \prod_{k=1}^n {\pi}_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta_{kT/n}-\theta_{(k-1)T/n}), \ \ \prod_{k=1}^n {\pi}_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta_{kT/n}-\theta_{(k-1)T/n}) \Big]}{Var \Big[ \prod_{k=1}^n {\pi}_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta_{kT/n}-\theta_{(k-1)T/n}) \Big]}. \end{align}$$(3.14) We note that $$\begin{align} E \left[ \prod_{k=1}^n {\pi}_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta_{kT/n}-\theta_{(k-1)T/n}) \right] =1, \end{align}$$(3.15) which is shown as follows: by the tower property of conditional expectation, we have $$\begin{equation} E \left[ \prod_{k=1}^n {\pi}_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta_{kT/n}-\theta_{(k-1)T/n}) \right] \end{equation}$$(3.16) $$\begin{align} = & E \left[ \prod_{k=1}^{n-1} {\pi}_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta_{kT/n}-\theta_{(k-1)T/n}) \Big. E \Big[{\pi}_{T/n}^{\bar{X}_{(n-1)T/n}^{(n),x}}(\theta_{T}-\theta_{(n-1)T/n}) \Big| \mathcal{F}_{(n-1)T/n} \Big] \right] \end{align}$$(3.17) $$\begin{equation} = E \left[ \prod_{k=1}^{n-1} {\pi}_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta_{kT/n}-\theta_{(k-1)T/n}) \right], \end{equation}$$(3.18) since it holds that $$\begin{align} \Big. E \Big[{\pi}_{T/n}^{\bar{X}_{(n-1)T/n}^{(n),x}}(\theta_{T}-\theta_{(n-1)T/n}) \Big| \mathcal{F}_{(n-1)T/n} \Big] =1 \ \ \text{a.s.} \end{align}$$(3.19) due to the representation of the weight given in Proposition 2.12 and the formula (3.9). Using the tower property of the conditional expectation repeatedly, we obtain $$\begin{align}E \left[ \prod_{k=1}^n {\pi}_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta_{kT/n}-\theta_{(k-1)T/n}) \right] &= E \left[ \prod_{k=1}^{n-1} {\pi}_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta_{kT/n}-\theta_{(k-1)T/n}) \right]\\ &\vdots \nonumber \\&= E \Big[{\pi}_{T/n}^{x}(\theta_{T/n}) \Big] =1. \end{align}$$(3.21) Then, the equation (3.15) results. We compute |$\lambda $| numerically from sample paths in the implementation since it does not usually have a closed form. Then, it is obvious that we can easily apply the Monte Carlo method to approximate the expectation |$E[Y_{(n)}^{CV1}]$| with a slight computational cost of computing |$\lambda $|⁠. Hence, instead of approximating the numerical scheme (3.1), it is more effective to compute the expectation |$E[Y_{(n)}^{CV1}]$| if the random variables |$f( \bar{X}_T^{(n),x} ) \prod _{k=1}^n {\pi }_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta _{kT/n}-\theta _{(k-1)T/n})$| and |$\prod _{k=1}^n {\pi }_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta _{kT/n}-\theta _{(k-1)T/n})$| are highly correlated. 3.2.2 Control variate 2 Next, we propose the second control variate as follows: $$\begin{align} Y_{(n)}^{CV2}: = f( \bar{X}_T^{(n),x} ) \prod_{k=1}^n {\pi}_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta_{kT/n}-\theta_{(k-1)T/n}) - \gamma \{ f( \bar{X}_T^{(n),x} ) - E [\,f( \bar{X}_T^{(n),x} ) ] \}, \end{align}$$(3.22) where |$\gamma \in \mathbf{R}$| is given by $$\begin{align} \gamma = \frac{Cov \Big[\,f( \bar{X}_T^{(n),x} ) \prod_{k=1}^n {\pi}_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta_{kT/n}-\theta_{(k-1)T/n}), \ \ f( \bar{X}_T^{(n),x} ) \Big]}{Var \Big[\,f( \bar{X}_T^{(n),x} ) \Big]}. \end{align}$$(3.23) Since in general it is hard to get an analytical formula for |$\gamma $|⁠, we derive numerically the value from the sample paths in the implementation. While it is expected that |$f( \bar{X}_T^{(n),x} )$| is highly correlated with |${f( \bar{X}_T^{(n),x} ) \prod _{k=1}^n {\pi }_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta _{kT/n}-\theta _{(k-1)T/n})}$|⁠, we have to know the exact value of |$E [\,f( \bar{X}_T^{(n),x} ) ]$| in order to apply the Monte Carlo method to |$E[ Y_{(n)}^{CV2}]$|⁠. Then, one may apply a quasi-Monte Carlo method to approximate |$E [\,f( \bar{X}_T^{(n),x} )]$| in the implementation of numerical simulation. Although there is an additional computational cost, the approximation of |$E [\,f( \bar{X}_T^{(n),x} )]$| with a quasi-Monte Carlo computing |$E[ Y_{(n)}^{CV2}]$|⁠, instead of |$(Q_{T/n})^n f(x) $|⁠, will lead to a variance reduction because of the correlation between $$\begin{equation*}f( \bar{X}_T^{(n),x} ) \ \ \mbox{and} \ \ {{f( \bar{X}_T^{(n),x} ) \prod_{k=1}^n {\pi}_{T/n}^{\bar{X}_{(k-1)T/n}^{(n),x}}(\theta_{kT/n}-\theta_{(k-1)T/n})}}.\end{equation*}$$ 4. Numerical results This section shows numerical results for the proposed scheme for various stochastic models. In Sections 4.1 and 4.2 a quasi-Monte Carlo method is applied to the proposed second-order scheme, and Section 4.3 gives numerical comparisons between our new scheme and another second-order scheme in the literature. Finally, Section 4.4 provides some results of standard Monte Carlo experiments for the new second-order scheme and the control variate methods introduced in the previous section. All numerical experiments are run on a Macbook Pro with a 2.3 GHz Intel Core i5 processor and 8 GB memory. 4.1 Asian option under SABR model: three-dimensional SDE First, we consider the following three-dimensional SDE for Asian option pricing under a stochastic volatility model: $$\begin{equation} {\hskip-100pt}d X_t^{x,1} = \; (X_t^{x,1})^{\beta} X_t^{x,2} \textrm{d}W_t^1, \ \ X_0^{x,1} = x_1>0, \end{equation}$$(4.1) $$\begin{equation} {\hskip-44pt}d X_t^{x,2} = \; \nu X_t^{x,2} ( \rho \textrm{d}W_t^1 + \sqrt{ 1- \rho^2 } \textrm{d}W_t^2), \ \ X_0^{x,2} = x_2>0, \end{equation}$$(4.2) $$\begin{equation} {\hskip-140pt}d X_t^{x,3} = \; X_t^{x,1} \textrm{d}t, \ \ X_0^{x,3} =x_3 \geqslant 0. \end{equation}$$(4.3) where |$\nu \geqslant 0,\; \beta \in [0,1]$| and |$\rho \in (-1,1)$|⁠. The process |$\{(X_t^{x,1},X_t^{x,2})\}_t$| is known as SABR (Stochastic Alpha Beta Rho) model, which is one of the important classes of stochastic volatility models in financial mathematics or/and in practice. Here, |$X_t^{x,1}$| and |$X_t^{x,2}$| represent the underlying asset and its volatility. In Asian option pricing, the time average of |$X^{x,1}$|⁠, i.e. |${\frac{1}{T}X_T^{x,3}=\frac{1}{T}\int _0^T X_t^{x,1} \textrm{d}t }$|⁠, is monitored. Our goal is to approximate the price of the Asian call option $$\begin{align} E[\,f_K ( X_T^x )], \end{align}$$(4.4) where |$ f_K ( \xi ) = \max ( \xi _3/T - K, 0),\; \xi \in \mathbf{R}^3$| with a strike price |$K>0$|⁠. We numerically compute |$E[\,f_K ( X_T^x )]$| by a quasi-Monte Carlo (in short QMC) method for various time steps. Then, we analyze the numerical error |$| E [\,f (X_T^x)] - (Q_{T/n})^n f(x)| $| and compare it with the error derived from the EM scheme. Here, we do not know the analytical solution for |$E[\,f_K ( X_T^x )], $| and hence compute the benchmark value by EM scheme with time steps |$n = 2^{11}$| and trajectories |$10^7$|⁠. We price the Asian call option by both approximation schemes with the parameters |$x_1 =100,\; x_2 = 0.3, \; x_3=0, \; \nu =0.1, \; \rho =-0.5, \; T=1$| in the cases |$\beta =0.5$| or |$0.3$|⁠. Under |$\beta =0.5$| we compute the expectation by QMC with trajectories |$10^6$|⁠, while we adopt |$2 \times 10^6$| paths under |$\beta =0.3$|⁠. Figure 1 (a) and (b) indicate the numerical errors |${\cal E}_{\mathrm{EM}}(n,K):=E [\,f_K( X_T^x ) ] - E[\,f_K( \bar{X}_T^{\mathrm{EM},(n),x})]$| for |$n=2^1,2^2,2^3$| and |${\cal E}_{\mathrm{Malliavin}}(n,K):=E [\,f_K( X_T^x) ] - (Q_{T/n})^n f_K(x)$| for |$n=2^0,2^1,2^2$| with strike prices |$K=10,\cdots ,200$|⁠. In particular, Figure 1 (a) and (b) show the numerical results with the parameter |$\beta =0.5$| and |$\beta =0.3$|⁠, respectively. In order to see the rate of convergence of both schemes, we check how the accuracy (in terms of maximum error with respect |$K$|⁠) will be improved as the number of time steps increases. Figure 2 shows the absolute errors |$|{\cal E}_{\mathrm{EM}}(n,K^{\ast })|$| and |$|{\cal E}_{\mathrm{Malliavin}}(n,K^{\ast \ast })|$|⁠, where |${K}^{\ast }= \textrm{argmax}_{K=10, \ldots , 200} |{\cal E}_{\mathrm{EM}}(2,K)|$| and |${K}^{\ast \ast } = \textrm{argmax}_{K=10, \ldots , 200} |{\cal E}_{\mathrm{Malliavin}}(1,K)|$| with |$\beta = 0.5$|⁠. We can see that the proposed scheme behaves as achieving |$O(n^{-2})$|⁠; on the other hand, EM scheme does as |$O(n^{-1}).$| Through this experiment it turns out that the new scheme has a superior efficiency compared to the EM scheme and achieves second-order weak convergence. Fig. 1. Open in new tabDownload slide Numerical errors for Asian option price (QMC). Fig. 1. Open in new tabDownload slide Numerical errors for Asian option price (QMC). Fig. 2. Open in new tabDownload slide Convergence rate for Asian option price (⁠|$\beta =0.5$|⁠) under SABR model (QMC). Fig. 2. Open in new tabDownload slide Convergence rate for Asian option price (⁠|$\beta =0.5$|⁠) under SABR model (QMC). 4.2 VWAP option under the SABR model: five-dimensional SDE Next, let us consider the VWAP call option pricing under the SABR stochastic volatility model, which has become important in practice nowadays. For instance, see Stace (2007) and Novikov et al. (2014). Let us consider the following five-dimensional SDE defined as follows: $$\begin{equation} {\hskip-145pt} \mathrm{d} X_t^{x,1} = \; (X_t^{x,1})^{\beta} X_t^{x,2} \textrm{d}W_t^1, \ \ X_0^{x,1} = x_1>0,\: \end{equation}$$(4.5) $$\begin{equation} {\hskip-80pt}\mathrm{d} X_t^{x,2} = \; \nu_1 X_t^{x,2} ( \rho_1 \textrm{d}W_t^1 + \sqrt{ 1- \rho_1^2} \textrm{d}W_t^2), \ \ X_0^{x,2} = x_2>0,\: \end{equation}$$(4.6) $$\begin{equation} \mathrm{d} X_t^{x,3} =\; (\theta -\kappa X_t^{x,3})\textrm{d}t + \nu_2 (X_t^{x,3})^{\lambda}( \rho_2 \textrm{d}W_t^1 + \sqrt{1- \rho_2^2}\textrm{d}W_t^3), \ \ X_0^{x,3} = x_3 \geqslant 0,\; \end{equation}$$(4.7) $$\begin{equation}{\hskip-169pt} \mathrm{d} X_t^{x,4} =\; X_t^{x,1} X_t^{x,3} \textrm{d}t, \ \ X_0^{x,4} = x_4 \geqslant0, \end{equation}$$(4.8) $$\begin{equation}{\hskip-183pt} \mathrm{d} X_t^{x,5} = \; X_t^{x,3} \textrm{d}t, \ \ X_0^{x,5} =x_5 \geqslant 0,\; \end{equation}$$(4.9) where |$\beta , \lambda \in [0,1],\; \nu _1, \nu _2 \geqslant 0, \; \rho _1, \rho _2 \in (-1,1)$| and |$\theta , \kappa $| are some constants. |$X_t^{x,1}$| and |$X_t^{x,2}$| represent an underlying asset price and that volatility respectively as in the previous example. |$X_t^{x,3}$| means the trading volume of stock and the trading volume weighted average of stock |$\int _0^T X_t^{x,1} X_t^{x,3} \textrm{d}t / \int _0^T X_t^{x,3} \textrm{d}t$| is the input of payoff function. Here we aim to approximate the option price |$E[\,f_K (X_T^x)] $| where |$f_K ( \xi ) = \max ( {\xi ^4}/{\xi ^5} - K, 0),\; \xi \in \mathbf{R}^5$| with a strike price |$K>0$|⁠. Then, we discuss the error between |$E[\,f_K (X_T^x)]$| and the proposed scheme. We set |$x_1=100,\;x_2=0.3, \;x_3=0,\; x_4= 0,\; x_5 = 0, \beta = \lambda =0.5, \; \rho_1 = -0.5, \; \rho _2 = 0.3, \; \nu _1= 0.1,\; \nu _2=0.3,\; \theta =10, \; \kappa =0.1, \; T=1.$| We derive the benchmark value by EM scheme with time steps |$n = 2^{11}$| and sample paths |$10^7$|⁠. Figure 3 (a) shows the numerical errors |${\cal E}_{\mathrm{EM}}(n,K):=E [\,f_K( X_T^x) ] - E[\,f_K( \bar{X}_T^{\mathrm{EM}, (n), x})]$| for |$n=2^1,2^2,2^3$| and |${\cal E}_{\mathrm{Malliavin}}(n,K):=E [\,f_K( X_T^x) ] - (Q_{T/n})^n f_K(x)$| for |$n=2^0,2^1,2^2$| with strike prices |$K=10,20, \ldots ,200$|⁠, where the approximation schemes are computed by QMC with sample paths |$10^6$|⁠. As is the case of Asian option under the SABR model, we can see that the new scheme outperforms the EM scheme uniformly even if the discretization step is small. In Figure 3 (b) we plot the absolute errors |$|{\cal E}_{\mathrm{EM}}(n,K^{\ast }) |$| and |$| {\cal E}_{\mathrm{Malliavin}}(n,K^{\ast \ast }) |$|⁠, where |${K}^{\ast }= \textrm{argmax}_{K=10, \ldots , 200}|{\cal E}_{\mathrm{EM}}(2,K) |$| and |$K^{\ast \ast } = \textrm{argmax}_{K=10, \ldots , 200} | {\cal E}_{\mathrm{Malliavin}}(1,K) | $|⁠. From this figure it is obvious that the proposed scheme achieves second-order weak approximation while EM scheme has first order. Fig. 3. Open in new tabDownload slide Numerical errors for VWAP option price (QMC). Fig. 3. Open in new tabDownload slide Numerical errors for VWAP option price (QMC). 4.3 Comparison of the new scheme and other second-order scheme In this subsection we compare our proposed scheme to another second-order weak approximation scheme in the literature. As an alternative method we adopt Ninomiya–Victoir (NV) method and show numerical results computed with the NV method and the proposed scheme. See the details of the NV scheme in Ninomiya & Victoir (2008). 4.3.1 Asian option under the CIR model: two-dimensional SDE As a first example let us consider the following two-dimensional SDE: $$\begin{equation} {} \mathrm{d} X_t^{x,1} = \; \alpha (\theta - X_t^{x,1}) \textrm{d}t + \sigma \sqrt{X_t^{x,1}} \textrm{d}W_t^1, \ \ X_0^{x,1} =x_1> 0, \end{equation}$$(4.10) $$\begin{equation} {}\mathrm{d} X_t^{x,2} = \; X_t^{x,1} \textrm{d}t, \ \ X_0^{x,2} = x_2 \geqslant 0, \end{equation}$$(4.11) where |$\alpha ,\theta $| and |$\sigma $| are some positive constants satisfying |$2 \alpha \theta -\sigma ^2>0$| to ensure the existence of the solution of SDE. This model is called Cox–Ingersoll–Ross (CIR) model in financial mathematics. Here, we also compute the price of the Asian call option |$E [\,f_K ( X_T^x )]$| where |$f_K ( \xi ) = \max ( \xi _2 / T - K), \; \xi \in \mathbf{R}^2$| with a strike price |$K>0$|⁠. Then, we numerically compute the expectation by using our proposed scheme and the NV method with QMC (⁠|$10^6$| paths). We set |$x_1= 100, \; x_2 =0, \; \alpha = 0.1, \; \ =100, \; \sigma = 2$| and |$T=1$| and derive the benchmark value by the EM scheme with time steps |$n= 2^{11}$| and paths |$10^7$|⁠. Figure 4 (a) shows the numerical errors |$\mathcal{E}_{\mathrm{Malliavin}}(n, K):= | E[\,f_K(X_T^x)] - (Q_{T/n})^n f(x)|$| and |$\mathcal{E}_{\mathrm{NV}}(n, K):= | E[\,f_K(X_T^x)] - (Q_{T/n}^{\mathrm{NV}})^n f(x)|$|⁠, where the operator |$\{ Q_t^{\mathrm{NV}} \}_{t>0} $| is Ninomiya–Victoir scheme. Figure 4 (b) plots the absolute errors |$|\mathcal{E}_{\mathrm{EM}}(n, K^{\ast })|, |\mathcal{E}_{\mathrm{Malliavin}}(n, K^{\ast \ast })|$| and |$|\mathcal{E}_{\mathrm{NV}}(n, K^{\ast \ast \ast })|$|⁠, where |${K}^{\ast }= \textrm{argmax}_{K=10, \ldots , 200}|{\cal E}_{\mathrm{EM}}(2,K) |, {K}^{\ast \ast }= \textrm{argmax}_{K=10, \ldots , 200}|{\cal E}_{\mathrm{Malliavin}}(1,K) |\ \mathrm{ and} \ \ K^{\ast \ast \ast } = \textrm{argmax}_{K=10, \ldots , 200} | {\cal E}_{\mathrm{NV}}(1,K) |$|⁠. From these figures we see that the both second-order schemes behave as |$O(n^{-2})$|⁠, while the EM scheme does as |$O (n^{-1})$|⁠. It is remarkable that the new second-order scheme provides a better approximation than the NV scheme for each time step |$n=1,2,4,8$| because of the better estimation when |$n=1$|⁠. Fig. 4. Open in new tabDownload slide Numerical errors for Asian option price under CIR model (QMC). Fig. 4. Open in new tabDownload slide Numerical errors for Asian option price under CIR model (QMC). Now, let us compare the computational times for the second-order schemes. Table 1 shows the accuracy and computational time of three weak approximation schemes (Euler–Maruyama, Ninomiya–Victoir and the new scheme). The label ‘Malliavin’ represents the new scheme. The parameters of the results are the same as those given in Figure 4 (b), except for the number of time steps and the number of paths. From Table 1 we observe that the computational time of the new scheme is less than half that of that for the NV scheme. This is because the NV scheme requires a |$(d+2)$|-fold composition of solutions of ODEs in each one-step simulation, while the new scheme (Malliavin) can be implemented with a single computation of the Malliavin weight |$\pi _{T/n}^{\bar{X}^{(n)}((k-1)T/n,x)}(\theta _{kT/n}- \theta _{(k-1)T/n})$| and |$\bar{X}^{(n)}(kT/n,x)$| in the one-step simulation. The same fact is observed in other SDEs, for instance in Table 2 in Section 4.3.2. Table 1 The accuracy and the computational time of Asian option price under CIR model Method . # of time steps |$n$| . # of paths |$M$| . Absolute errors . Time (s) . Euler–Maruyama 256 400000 0.02717 1.27 Ninomiya–Victoir 4 200000 0.03259 0.57 Malliavin 4 200000 0.01498 0.27 Method . # of time steps |$n$| . # of paths |$M$| . Absolute errors . Time (s) . Euler–Maruyama 256 400000 0.02717 1.27 Ninomiya–Victoir 4 200000 0.03259 0.57 Malliavin 4 200000 0.01498 0.27 Open in new tab Table 1 The accuracy and the computational time of Asian option price under CIR model Method . # of time steps |$n$| . # of paths |$M$| . Absolute errors . Time (s) . Euler–Maruyama 256 400000 0.02717 1.27 Ninomiya–Victoir 4 200000 0.03259 0.57 Malliavin 4 200000 0.01498 0.27 Method . # of time steps |$n$| . # of paths |$M$| . Absolute errors . Time (s) . Euler–Maruyama 256 400000 0.02717 1.27 Ninomiya–Victoir 4 200000 0.03259 0.57 Malliavin 4 200000 0.01498 0.27 Open in new tab Table 2 The accuracy and the computational time of Asian option price under Heston model Method . # of time steps |$n$| . # of paths |$M$| . Absolute errors . Time (s) . Euler–Maruyama 256 400000 0.02303 4.87 Ninomiya–Victoir 4 200000 0.02727 2.43 Malliavin 4 200000 0.00629 0.32 Method . # of time steps |$n$| . # of paths |$M$| . Absolute errors . Time (s) . Euler–Maruyama 256 400000 0.02303 4.87 Ninomiya–Victoir 4 200000 0.02727 2.43 Malliavin 4 200000 0.00629 0.32 Open in new tab Table 2 The accuracy and the computational time of Asian option price under Heston model Method . # of time steps |$n$| . # of paths |$M$| . Absolute errors . Time (s) . Euler–Maruyama 256 400000 0.02303 4.87 Ninomiya–Victoir 4 200000 0.02727 2.43 Malliavin 4 200000 0.00629 0.32 Method . # of time steps |$n$| . # of paths |$M$| . Absolute errors . Time (s) . Euler–Maruyama 256 400000 0.02303 4.87 Ninomiya–Victoir 4 200000 0.02727 2.43 Malliavin 4 200000 0.00629 0.32 Open in new tab In addition, we analyze the variance of both second-order schemes, applying the standard Monte-Carlo (MC in short) method. Here, we approximate |$E[\,f_K (X_T^x)]$| by the NV scheme and our proposed scheme using MC with paths |$5 \times 10^4$|⁠, and repeat the computations |$1000$| times to derive the standard deviation. The parameters are given as follows: |$x_1= 100, \; x_2 =0, \; \alpha = 0.1, \; \theta =100, \; \sigma = 2, \; K=100$| and |$T=1$|⁠. Furthermore, we use the control variate methods introduced in Section 3.2 (control variate 1 and control variate 2) to see the effect of variance reduction. In Figure 5, as |$y$|-axis, we plot the absolute error between the benchmark value and the average of |$1000$| results of MC computation, while the standard deviation as the |$x$|-axis. The labels (Malliavin + control variate 1) and (Malliavin + control variate 2) mean |$E[ Y_{(n)}^{CV1}]$| (in (3.13)) and |$E[ Y_{(n)}^{CV2}]$| (in (3.22)). Also, the label |$(n),\; n=1,2,4,8$|⁠, indicates the number of time steps. From this figure we see that when |$n=1$|⁠, the NV scheme gives a smaller variance than the proposed scheme; however, both schemes achieve almost the same variance as the number of time steps reaches |$8$|⁠. We also notice from Fig. 5 that the control variate methods work for our proposed scheme and, especially, the second control variate method reduces dramatically the standard deviation compared to the other three approximation schemes. Fig. 5. Open in new tabDownload slide Standard deviation of Asian option pricing under CIR model (MC). Fig. 5. Open in new tabDownload slide Standard deviation of Asian option pricing under CIR model (MC). 4.3.2 Asian option under Heston model: three-dimensional SDE Next, we consider the following three-dimensional SDE called Heston model in financial mathematics: $$\begin{equation} \mathrm{d} X_t^{x,1} =\; X_t^{x,1} \sqrt{X_t^{x,2}} \textrm{d}W_t^1, \ \ X_0^{x,1} = x_1>0, \end{equation}$$(4.12) $$\begin{equation} \mathrm{d} X_t^{x,2} = \; \alpha (\theta - X_t^{x,2}) \textrm{d}t + \beta \sqrt{X_t^{x,2}} (\rho \textrm{d}W_t^1 \ + \sqrt{1- \rho^2} \textrm{d}W_t^2), \ \ X_0^{x,2} = x_2>0, \end{equation}$$(4.13) $$\begin{equation} \mathrm{d} X_t^{x,3} = \; X_t^{x,1} \textrm{d}t, \ \ X_0^{x, 2} = x_3 \geqslant 0, \end{equation}$$(4.14) where |$\alpha ,\theta $| and |$\beta $| are some positive constants and |$\rho \in (-1,1)$|⁠. Then, we set the parameters as follows: |$x_1 =100, \; x_2 =0.09, \; x_3 =0, \; \alpha = 2.0, \theta = 0.09, \; \beta = 0.1, \; \rho =0.7, \; T=1$|⁠. Here, again we consider the Asian call option price |$E[\,f_K ( X_T^x)],$| where |$f_K(\xi ) = \max ( \xi _3 / T -K, 0)$| for |$\xi \in \mathbf{R}^3$| with a strike price |$K>0$|⁠. As we did in the previous subsection (CIR model), we approximately compute |$E[\,f_K ( X_T^x)]$| with NV scheme and our proposed scheme by using QMC (⁠|$2\times 10^6$| paths). The benchmark value of |$E[\,f_K ( X_T^x)]$| is derived from the Euler–Maruyama scheme with number of time steps |$n=2^{11}$| and paths |$10^7$|⁠. In Figure 6 (a) we plot the numerical errors |$\mathcal{E}_{\mathrm{Malliavin}}(n, K):= | E[\,f_K(X_T^x)] - (Q_{T/n})^n f(x)|$| and |$\mathcal{E}_{\mathrm{NV}}(n, K):= | E[\,f_K(X_T^x)] - (Q_{T/n}^{\mathrm{NV}})^n f(x)|$| for |$n=1,2,4$| and |$K=10, 20, $||$\ldots , 200$|⁠. On the other hand, Figure 6 (b) presents the absolute errors |$|\mathcal{E}_{\mathrm{EM}}(n, K^{\ast })|, |\mathcal{E}_{\mathrm{Malliavin}}(n, K^{\ast \ast })|$| and |$|\mathcal{E}_{\mathrm{NV}}(n, K^{\ast \ast \ast })|$|⁠, where |$K^{\ast }, K^{\ast \ast }$| and |$K^{\ast \ast \ast }$| are defined in Subsection 4.3.1. Similarly to the example of Asian option price under the CIR model, the maximum error (with respect to strike prices) derived from the proposed scheme is smaller than that of the NV scheme for every time step |$n=1,2, 4, 8$|⁠. Fig. 6. Open in new tabDownload slide Numerical errors for Asian option price under Heston model (QMC). Fig. 6. Open in new tabDownload slide Numerical errors for Asian option price under Heston model (QMC). Furthermore, we will mention the computational times for both of the second-order schemes. Table 2 shows the accuracy and computational time of the three weak approximation schemes (Euler–Maruyama, Ninomiya–Victoir and the new scheme). The parameters of the results are the same as the ones given in Fig. 6 (b), except for the number of time steps and the number of paths. The label ‘Malliavin’ represents the new scheme. From Table 2 we see that the NV scheme and our proposed scheme achieve the desiered accuracies with a fewer number of discretizations (⁠|$n=2$|⁠) compared to EM the scheme, and that the new scheme computes the price much faster than the NV scheme. 4.4 Variance analysis of the scheme and the effect of control variate In this subsection we verify our discussion in Section 3 through numerical experiments using the standard Monte Carlo method. Here, we consider two examples: mean energy of synchrotron oscillations and Asian call option pricing under the SABR model. In the first example we consider a degenerate and interacting system of two-dimensional SDE given by $$\begin{equation} \textrm{d}X_t^{x,1} = -w^2 \sin(X_t^{x,2})\textrm{d}t - \sigma_1 \cos(X_t^{x,2}) \textrm{d}W_t^1 - \sigma_2 \sin(X_t^{x,2}) \textrm{d}W_t^2, \ \ X_0^{x,1} =x_1, \end{equation}$$(4.15) $$\begin{equation} \textrm{d}X_t^{x,2} = X_t^{x,1} \textrm{d}t, \ \ X_0^{x,2} =x_2, \end{equation}$$(4.16) and the expectation |$E [ \mathcal{E}( X_T^x )]$| where |$\mathcal{E} (x) = (x_1)^2/2 - w^2 \cos (x_2), x \in \mathbf{R}^2$|⁠. This expectation is interpreted as the mean energy of the system. See the details in Seesselbeg et al. (1994) and Milstein & Tretyakov (2004). In the numerical experiments of |$E [ \mathcal{E}( X_T^x )]$| we set the parameters of the SDE as |$x_1 = 1, x_2 =0, w =2, \sigma _1 = \sigma _2 =0.2$| and |$T=1$|⁠. For the case of Asian option pricing under the SABR model, we adopt the same parameter setting given in Section 4.1. First, we apply the EM scheme and the proposed second-order scheme to approximate the expectation |$E [ \mathcal{E}( X_T^x )]$|⁠. Secondly, we also use those schemes to approximate (4.4) with |$K=K^{*}$| for EM scheme and |$K=K^{**}$| for the new scheme, where |$K^*$| and |$K^{**}$| are the same ones defined in Section 4.1. In order to analyze the variance of our new scheme for each time step, we calculate the standard deviation based on |$1000$| batches of a MC simulation with trajectories |$5 \times 10^4$| for |$E [ \mathcal{E}( X_T^x )]$| and |$10^5$| for (4.4). Under the same situation we also apply the MC method to the control variate methods (control variate 1 and control variate 2 introduced in Section 3.2) for the proposed second-order scheme in order to see the effect of variance reduction. In Figure 7 (a) and (b) we plot pairs of the absolute error and standard deviation for the expectations |$E [ \mathcal{E}( X_T^x )]$| and (4.4), respectively. The labels (EM), (Malliavin), (Malliavin + control variate 1) and (Malliavin + control variate 2) mean the numerical results by the EM scheme, the proposed second-order scheme, |$E[ Y_{(n)}^{CV1}]$| (in (3.13)) and |$E[ Y_{(n)}^{CV2}]$| (in (3.22)). The subscript |$(n), n=1,2, \ldots , 128$|⁠, means the number of discretizations. From these figures we find that in both cases the standard deviations of the new second-order scheme are basically stable as |$n$| increases. It is remarkable that in the Asian call option pricing the standard deviation of the new scheme does not rise, although we observe the increase of variance by the EM scheme as |$n$| rises. Moreover, in both graphs, the control variates succeed to reduce the variance of the numerical scheme. In particular, in the example of synchrotron oscillation, the standard deviations have been significantly reduced by the control variate |$Y_{(n)}^{CV1}$|⁠, while in the case of Asian call option pricing, the control variate |$Y_{(n)}^{CV2}$| leads to an effective variance reduction and the standard deviations are much lower than those of EM the scheme. Fig. 7. Open in new tabDownload slide Standard deviation of numerical schemes. Fig. 7. Open in new tabDownload slide Standard deviation of numerical schemes. Table 3 (a) and (b) illustrate mean squared errors (MSEs) derived from the EM scheme, the new scheme and two types of control variate. For both examples, the MSEs tend to be stable as the number of time steps |$n$| increases, and the proposed scheme gives a lower MSEs than those by the EM scheme. Moreover, due to the reduction of variance, the MSEs take lower values when we use the control variates. Table 3 Mean squared errors n . EM . Malliavin . Malliavin + control variate 1 . Malliavin + control variate 2 . 2 1.8091661247 0.3212905669 0.0165676026 0.3251797242 4 0.4378102970 0.0462449638 0.0004075337 0.0455611190 8 0.0865344672 0.0019798735 0.0000061083 0.0018201738 16 0.0181612771 0.0002238882 0.0000005630 0.0001833357 32 0.0041106952 0.0000470028 0.0000002999 0.0000362555 (a) energy of synchrotron oscillation (MC) 1 - 0.0290062249 0.0283413459 0.0257984965 2 0.6193434455 0.0049031599 0.0047323448 0.0023213852 4 0.2241217967 0.0029694654 0.0027570018 0.0006388366 8 0.0654723371 0.0027890036 0.0026865940 0.0003683108 16 0.0186115832 0.0026769685 0.0025203560 0.0004097712 (b) Asian call option pricing (MC) n . EM . Malliavin . Malliavin + control variate 1 . Malliavin + control variate 2 . 2 1.8091661247 0.3212905669 0.0165676026 0.3251797242 4 0.4378102970 0.0462449638 0.0004075337 0.0455611190 8 0.0865344672 0.0019798735 0.0000061083 0.0018201738 16 0.0181612771 0.0002238882 0.0000005630 0.0001833357 32 0.0041106952 0.0000470028 0.0000002999 0.0000362555 (a) energy of synchrotron oscillation (MC) 1 - 0.0290062249 0.0283413459 0.0257984965 2 0.6193434455 0.0049031599 0.0047323448 0.0023213852 4 0.2241217967 0.0029694654 0.0027570018 0.0006388366 8 0.0654723371 0.0027890036 0.0026865940 0.0003683108 16 0.0186115832 0.0026769685 0.0025203560 0.0004097712 (b) Asian call option pricing (MC) Open in new tab Table 3 Mean squared errors n . EM . Malliavin . Malliavin + control variate 1 . Malliavin + control variate 2 . 2 1.8091661247 0.3212905669 0.0165676026 0.3251797242 4 0.4378102970 0.0462449638 0.0004075337 0.0455611190 8 0.0865344672 0.0019798735 0.0000061083 0.0018201738 16 0.0181612771 0.0002238882 0.0000005630 0.0001833357 32 0.0041106952 0.0000470028 0.0000002999 0.0000362555 (a) energy of synchrotron oscillation (MC) 1 - 0.0290062249 0.0283413459 0.0257984965 2 0.6193434455 0.0049031599 0.0047323448 0.0023213852 4 0.2241217967 0.0029694654 0.0027570018 0.0006388366 8 0.0654723371 0.0027890036 0.0026865940 0.0003683108 16 0.0186115832 0.0026769685 0.0025203560 0.0004097712 (b) Asian call option pricing (MC) n . EM . Malliavin . Malliavin + control variate 1 . Malliavin + control variate 2 . 2 1.8091661247 0.3212905669 0.0165676026 0.3251797242 4 0.4378102970 0.0462449638 0.0004075337 0.0455611190 8 0.0865344672 0.0019798735 0.0000061083 0.0018201738 16 0.0181612771 0.0002238882 0.0000005630 0.0001833357 32 0.0041106952 0.0000470028 0.0000002999 0.0000362555 (a) energy of synchrotron oscillation (MC) 1 - 0.0290062249 0.0283413459 0.0257984965 2 0.6193434455 0.0049031599 0.0047323448 0.0023213852 4 0.2241217967 0.0029694654 0.0027570018 0.0006388366 8 0.0654723371 0.0027890036 0.0026865940 0.0003683108 16 0.0186115832 0.0026769685 0.0025203560 0.0004097712 (b) Asian call option pricing (MC) Open in new tab Remark 4.1 In Section 3.1 we gave the upper bound for the |$L^2$|-norm of the proposed second-order scheme to estimate the standard deviation. We especially confirm by (3.2) that as |$n$| becomes large the variance will never explode since the |$L^2$|-norm has an exponential bound independent of |$n$|⁠. However, the numerical experiments shown above imply that the variance of the new scheme tends to decrease as the number of time steps increases, and hence the bound (3.2) might still be crude. The derivation of an appropriate estimate for that is our future work. Remark 4.2 The SABR ((involving CEV, CIR)) type models do not satisfy the assumption (smoothness) on |$V_i$| in Section 2 rigorously since the diffusion coefficient has an unbounded derivative at zero. Despite this problem our scheme can be formally applicable and can work efficiently, which has been confirmed by the numerical examples. A rigorous treatment could be made using a smooth modification technique in Takahashi (2015) as follows. The original SABR model is $$\begin{align*} \textrm{d}X_t^{x,1}&=(X_t^{x,1})^{\beta} X_t^{x,2} \textrm{d}W_t^{1}, \ \ X_0^{x,1}=x_1>0,\\ \textrm{d}X_t^{x,2}&=\nu X_t^{x,2} (\rho \textrm{d}W_t^{1}+\sqrt{1-\rho^2} \textrm{d}W^2_t), \ \ X_0^{x,2}=x_2=\sigma=\epsilon >0. \end{align*}$$ The modified SABR model is $$\begin{align*} \mathrm{d}\widetilde{X_t}^{x,1}&=\epsilon g(\widetilde{X_t}^{x,1}) Y_t \textrm{d}W_t^{1}, \ \ \widetilde{X_0}^{x,1}=x_1>0,\\ \mathrm{d}Y_t \ &=\nu Y_t (\rho \textrm{d}W_t^{1}+\sqrt{1-\rho^2 }\textrm{d}W^2_t), \ \ Y_0=1. \end{align*}$$ Here, the initial volatility |$\sigma $| is regarded as a small parameter |$\sigma =\epsilon \in (0,1)$|⁠, and |$g$| is a smooth modification function given by $$\begin{align*} g(\xi)=\frac{\psi(\xi-\epsilon_2)}{\psi(\xi-\epsilon_2)+\psi(\epsilon_1-\xi)} \xi^\beta, \ \ \ 0<\epsilon_2<\epsilon_1, \end{align*}$$ for some small |$\epsilon _1>0$|⁠, where |$\psi (\xi )=e^{-1/\xi }$| (⁠|$\xi>0$|⁠), |$\psi (\xi )=0$| (⁠|$\xi \leqslant 0$|⁠). Then, for a measurable (or a smooth) function |$f$| we can show $$\begin{align*} E \left[ \left| f(X_t^{x,1})-f \big(\widetilde{X_t}^{x,1}\big) \right| \right]=O(\epsilon^k), \ \ \ \mbox{for all} \ \ k\in \mathbf{N}, \end{align*}$$ by applying a large deviation type inequality. Therefore, the difference between |$f(X_t^{x,1})$| and |$f \big (\widetilde{X_t}^{x,1}\big )$| is negligible in an asymptotic sense. The same technique can be also applied to the Heston-type model. 5. Concluding remarks We introduced a second-order weak approximation of a multidimensional degenerate system of SDEs satisfying a certain Hörmander condition. The new scheme contains the Markov chain scheme, which is nondegenerate in Malliavin’s sense and a Malliavin weight, which is explicitly given as the polynomial of Gaussian random variables (Section 2.5). Hence, we are able to apply the Monte Carlo method to the scheme for any SDE satisfying the assumptions in this paper. Furthermore, we showed with the help of Malliavin calculus that the new second-order scheme works with nonsmooth test functions. The numerical experiments for hypoelliptic diffusions suggested the validity of the scheme and the stability of the variance of scheme as the number of time steps increases. It is worth considering how the second-order scheme works under a weaker mathematical condition. For example, could the result be extended to the case when |$f$| is measurable and has a polynomial growth as in Section 6 in Bally & Talay (1996) using Watanabe distribution theory in our case. Also, the boundedness of the coefficients of SDEs might be relaxed. In particular, it is important to construct a second-order scheme under a weaker Hörmander condition. Application to nonlinear models such as backward SDEs may be also interesting. These should be our future tasks. Acknowledgements We are grateful to the editor, the associate editor and two anonymous reviewers for their valuable comments and suggestions. Funding JSPS KAKENHI (16K13773, 19K13736); MEXT, Japan; research fund from Tokio Marine Kagami Memorial Foundation. References Bally , V. ( 2003 ) An elementary introduction to Malliavin calculus . INRIA . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Bally , V. & Kohatsu-Higa , A. ( 2015 ) A probabilistic interpretation of the parametrix method . Ann. Appl. Probab. , 25 , 3095 – 3138 . Google Scholar Crossref Search ADS WorldCat Bally , V. & Talay , D. ( 1996 ) The law of the Euler scheme for stochastic differential equations I. Convergence rate of the distribution function . Probab. Theory Related Fields , 104 , 43 – 60 . Google Scholar Crossref Search ADS WorldCat Bayer , C. , Friz , P. & Loeffen , R. ( 2013 ) Semi-closed form cubature and applications to financial diffusion models . Quant. Finance , 13 , 769 – 782 . Google Scholar Crossref Search ADS WorldCat Beskos , A. & Roberts , G. O. ( 2005 ) Exact simulation of diffusions . Ann. Appl. Probab. , 15 , 2422 – 2444 . Google Scholar Crossref Search ADS WorldCat Gyurkó , L. G. & Lyons , T. ( 2010a ) Efficient and practical implementations of cubature on Wiener space . Stochastic Analysis 2010 . Springer , pp. 73 – 111 . Berlin . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Gyurkó , L. G. & Lyons , T. ( 2010b ) Rough paths based numerical algorithms in computational finance . Contemporary Mathematics , 515 . AMS , pp. 17 – 46 . Rhode Island . Google Scholar Crossref Search ADS WorldCat Hairer , M. , Stuart , A. M. & Voss , J. ( 2011 ) Sampling conditioned hypoelliptic diffusions . Ann. Appl. Probab. , 21 , 669 – 698 . Google Scholar Crossref Search ADS WorldCat Ikeda , N. & Watanabe , S. ( 1989 ) Stochastic Differential Equations and Diffusion Processes , 2nd edn. North-Holland Mathematical Library . Tokyo . Kloeden , P. E. & Platen , E. ( 1999 ) Numerical Solution of Stochastic Differential Equations . Springer . Berlin . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Kolmogorov , A. ( 1934 ) Zufällige Bewegungen . Ann. Math. (2) , 35 , 116 – 117 . Google Scholar Crossref Search ADS WorldCat Kunitomo , N. & Takahashi , A. ( 1992 ) Pricing average options (in Japanese) . Japan Financ. Rev. , 14 , 1 – 20 . Google Scholar OpenURL Placeholder Text WorldCat Kusuoka , S. ( 2001 ) Approximation of expectation of diffusion process and mathematical finance . Adv. Stud. Pure Math. , 31 , 147 – 165 . Google Scholar OpenURL Placeholder Text WorldCat Kusuoka , S. ( 2004 ) Approximation of expectation of diffusion process based on lie algebra and Malliavin calculus . Adv. Math. Econ. , 6 , 69 – 83 . Google Scholar Crossref Search ADS WorldCat Kusuoka , S. & Stroock , D. ( 1985 ) Applications of the Malliavin calculus part II . J. Fac. Sci. Univ. Tokyo , 32 , 1 – 76 . Google Scholar OpenURL Placeholder Text WorldCat Litterer , C. & Lyons , T. ( 2012 ) High order recombination and an application to cubature on wiener space . Ann. Appl. Probab. , 22 , 1301 – 1327 . Google Scholar Crossref Search ADS WorldCat Lyons , T. & Victoir , N. ( 2004 ) Cubature on Wiener space . Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. , 460 , 169 – 198 . Google Scholar Crossref Search ADS WorldCat Milstein , G. & Tretyakov , M. ( 2004 ) Stochastic Numerics for Mathematical Physics . Springer . Berlin . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Ninomiya , S. & Victoir , N. ( 2008 ) Weak approximation of stochastic differential equations and application to derivative pricing . Appl. Math. Finance , 15 , 107 – 121 . Google Scholar Crossref Search ADS WorldCat Novikov , A. A. , Ling , T. G. & Kordzakhia , N. ( 2014 ) Pricing of volume-weighted average options: analytical approximations and numerical results . Inspired by Finance (Y. Kabanov, M. Rutkowski & T. Zariphopoulou eds). Springer , 461 – 474 . Switzerland . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Nualart , D. ( 2006 ) The Malliavin Calculus and Related Topics . Springer . Berlin Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Pagès , G. ( 2018 ) Numerical Probability . Springer . Switzerland . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Seeßelberg , M. , Breuer , H. P., Mais , H., Petruccione , F. & Honerkamp , J. ( 1994 ) Simulation of one-dimensional noisy Hamiltonian systems and their application to particle storage rings . Z. Phys. C , 62 , 62 – 73 . Google Scholar Crossref Search ADS WorldCat Shigekawa , I. ( 2004 ) Stochastic Analysis . American Mathematical Society . Rhode Island . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Stace , A. W. ( 2007 ) A moment matching approach to the valuation of a volume weighted average price option . Int. J. Theor. Appl. Finance , 10 , 95 – 110 . Google Scholar Crossref Search ADS WorldCat Takahashi , A. ( 2015 ) Asymptotic expansion approach in finance . Large Deviations and Asymptotic Methods in Finance (P. Friz, J. Gatheral, A. Gulisashvili, A. Jacquier & J. Teichmann eds). Springer Proceedings in Mathematics & Statistics . Switzerland . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Takahashi , A. & Yamada , T. ( 2012 ) An asymptotic expansion with push-down of Malliavin weights . SIAM J. Financ. Math. , 3 , 95 – 136 . Google Scholar Crossref Search ADS WorldCat Takahashi , A. & Yamada , T. ( 2016 ) A weak approximation with asymptotic expansion and multidimensional Malliavin weights . Ann. Appl. Probab. , 26 , 818 – 856 . Google Scholar Crossref Search ADS WorldCat Watanabe , S. ( 1987 ) Analysis of wiener functionals (Malliavin calculus) and its applications to heat kernels . Ann. Probab. , 15 , 1 – 39 . Google Scholar Crossref Search ADS WorldCat Yamada , T. ( 2019 ) An arbitrary high order weak approximation of SDE and Malliavin Monte Carlo: analysis of probability distribution functions . SIAM J. Numer. Anal. , 57 , 563 – 591 . Google Scholar Crossref Search ADS WorldCat Yamada , T. & Yamamoto , K. ( 2018a ) A second-order discretization with Malliavin weight and Quasi-Monte Carlo simulation . Quant. Finance (Published Online) . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Yamada , T. & Yamamoto , K. ( 2018b ) A second order weak approximation of SDEs using Markov chain without Lévy area simulation . Monte Carlo Methods Appl. , 24 , 289 – 308 . Google Scholar Crossref Search ADS WorldCat Appendix A. Proof of Lemma 2.4 By the stochastic Taylor expansion, we have $$\begin{eqnarray*} X_t^{x,i_1,\varepsilon}&=&x_{i_1}+\varepsilon^2 V_0^{i_1}(x)t+\varepsilon \sum_{j=1}^d V_j^{i_1}(x)W^{j}_t\\ && \ \ +\sum_{k=2}^m \varepsilon^k \sum_{J=(j_1,\cdots,j_l),\|J\|=k,|J|\geqslant 2} \hat{V}_{j_{1}} \cdots \hat{V}_{j_{l-1}}V_{j_l}^{i_1}(x)I_{(j_1,\cdots,j_l)}(t)+\varepsilon^{m+1} R^{i_2}_{m+1}(\varepsilon), \\ X_t^{x,i_2,\varepsilon}&=&x_{i_2}+\varepsilon^2 V_0^{i_2}(x)t+\varepsilon^3 \sum_{i=1}^d \int_0^t \int_0^s \hat{V}_{j}V_0^{i_2}( X_u^{x,\varepsilon} ) \textrm{d}W_u^i \textrm{d}s+\varepsilon^4 \int_0^t \int_0^s \hat{V}_{0}V_0^{i_2}( X_u^{x,\varepsilon} ) \textrm{d}u \textrm{d}s\\ &=&x_{i_2}+\varepsilon^2 V_0^{i_2}(x)t+\varepsilon^3 \sum_{j=1}^d \hat{V}_{j}V_0^{i_2}(x)\int_0^t W^{j}_s \textrm{d}s + \varepsilon^4 \int_0^t \int_0^s \hat{V}_0 V_0^{i_1}(x) \textrm{d}u \textrm{d}s \\ && \ \ +\sum_{k=2}^m \varepsilon^{k+2} \sum_{J=(j_1,\cdots,j_l), \|J\|=k,|J|\geqslant2} \hat{V}_{j_{1}} \cdots \hat{V}_{j_{l}}V_{0}^{i_2}(x) I_{(j_1,\cdots,j_l,0)}(t)+\varepsilon^{m+3} R^{i_2}_{m+1}(\varepsilon). \end{eqnarray*}$$ We define $$\begin{align} Y_t^{i_1,\varepsilon}:=&\;\frac{X_t^{x,i_1,\varepsilon}-x_{i_1}-\varepsilon^2 V_0^{i_1}(x)t}{\varepsilon} \nonumber \\ =&\; \sum_{j=1}^d V_j^{i_1}(x)W^{j}_t+\sum_{k=2}^m \varepsilon^{k-1} \sum_{{J=(j_1,\cdots,j_l),\|J\|=k,|J|\geqslant 2}} \hat{V}_{j_{1}} \cdots \hat{V}_{j_{l-1}}V_{j_l}^{i_1}(x)I_{(j_1,\cdots,j_l)}(t)+\varepsilon^{m} R^{i_1}_{m+1}(\varepsilon), \end{align}$$(A.1) $$\begin{align} Y_t^{i_2,\varepsilon}:=&\; \frac{X_t^{x,i_2,\varepsilon}-x_{i_2}-\varepsilon^2 V_0^{i_2}(x)t -\varepsilon^4 \hat{V}_0 V_0^{i_2}(x) \frac{t^2}{2}} {\varepsilon^3} \nonumber \\ =&\! \sum_{j=1}^d \hat{V}_{j}V_0^{i_2}(x)\int_0^t W^{j}_s \textrm{d}s +\sum_{k=2}^m \varepsilon^{k-1} \sum_{{J=(j_1,\cdots,j_l), \|J\|=k,|J|\geqslant2}} \hat{V}_{j_{1}} \cdots \hat{V}_{j_{l}}V_{0}^{i_2}(x) I_{(j_1,\cdots,j_l,0)}(t)+\varepsilon^{m} R^{i_2}_{m+1}(\varepsilon). \end{align}$$(A.2) For |$\varphi \in C^{\infty }_b( \mathbf{R}^N)$|⁠, it holds $$\begin{align} E[ \varphi (X_1^{x, \varepsilon}) ]= \int_{\mathbf{R}^N} \varphi ( \tau (y)) {}_{-\infty} \langle \delta_y (Y_1^{\varepsilon}),1 \rangle_{\infty} \textrm{d}y. \end{align}$$(A.3) We expand |$E[\delta _y (Y_1^{\varepsilon })]={}_{-\infty } \langle \delta _y (Y_1^{\varepsilon }),1 \rangle _{\infty }$| to obtain the expansion of |$E[\varphi (X_1^{\varepsilon })]$|⁠. We first check the following condition on the Malliavin covariance matrix |$Y_1^{\varepsilon }$|⁠: for |$p<\infty $| there is a constant |$C>0$| such that $$\begin{align} \sup_{\varepsilon \in [0,1]} \| (\det \sigma^{Y_1^{\varepsilon}})^{-1} \|_p \leqslant C, \end{align}$$(A.4) for all |$x \in \mathbf{R}^N$|⁠. Using stochastic Taylor expansion we have the expansion of |$(J_s^{\varepsilon })^{-1} \varepsilon V_{\alpha } ( X_s^{\varepsilon } )$|⁠, |$\alpha \in \{ 1, \ldots , d \}$| as follows: $$\begin{align*} (J_s^{\varepsilon})^{-1} \varepsilon V_{\alpha} ( X_s^{\varepsilon} ) =&\; \varepsilon V_{\alpha}(x) + \varepsilon^3 [ \widetilde{V}_0, {V}_{\alpha} ](x) s + R^{\varepsilon}(s, x), \end{align*}$$ where $$\begin{align*} R^{\varepsilon}(s, x)=&\; \varepsilon^5 \int_0^s \int_0^u (J_v^{\varepsilon})^{-1} [ \widetilde{V}_0, [ \widetilde{V}_0, {V}_{\alpha}]] (X_v^{\varepsilon}) \textrm{d}v \textrm{d}u + \varepsilon^4 \sum_{k=1}^d \int_0^s \int_0^u (J_v^{\varepsilon})^{-1} [ {V}_k, [ \widetilde{V}_0, {V}_{\alpha}]] (X_v^{\varepsilon}) \circ \textrm{d}W_v^k \textrm{d}u \\ &+ \varepsilon^2 \sum_{k=1}^d \int_0^s (J_u^{\varepsilon})^{-1} [ {V}_k, {V}_\alpha] ( X_u^{\varepsilon}) \circ \textrm{d}W_u^{k}. \end{align*}$$ We have for |$i_1=1,\ldots ,K$| and |$i_2 = K+1, \ldots , N$|⁠, $$\begin{align*} \frac{1}{\varepsilon} [(J_s^{\varepsilon})^{-1} \varepsilon V_{\alpha} ( X_s^{\varepsilon} ) ]^{i_1} &= V_{\alpha}^{i_1}(x) + \varepsilon^2 [ \widetilde{V}_0, {V}_{\alpha} ]^{i_1}(x) s + \widetilde{R}^{\varepsilon,i_1}(s, x)= c_{1, \alpha}^{\varepsilon, i_1}(s, x) +c_{2, \alpha}^{\varepsilon, i_1}(s, x)+ \widetilde{R}^{\varepsilon,i_1}(s, x), \\ \frac{1}{\varepsilon^3} [(J_s^{\varepsilon})^{-1} \varepsilon V_{\alpha} ( X_s^{\varepsilon} ) ]^{i_2} &= [ \widetilde{V}_0, {V}_{\alpha} ]^{i_2}(x) s + \widetilde{R}^{\varepsilon,i_2}(s, x) = c_{2,\alpha}^{\varepsilon, i_2}(s, x) + \widetilde{R}^{\varepsilon,i_2}(s, x), \end{align*}$$ where |$N$|-dimensional vectors |$c_{1,\alpha }^{\varepsilon }(s,x) $| and |$c_{2,\alpha }^{\varepsilon }(s,x), \; \alpha \in \{1,\ldots , d\}$| are defined as $$\begin{align*} c_{1,\alpha}^{\varepsilon}(s,x) & = (V_\alpha^1(x), \ldots, V_{\alpha}^K(x),0 \ldots,0), \\ c_{2, \alpha}^{\varepsilon}(s,x) &= ( \varepsilon^2 [\widetilde{V}_0, V_{\alpha}]^1(x)s, \ldots, \varepsilon^2[ \widetilde{V}_0, V_{\alpha}]^K(x) s, [\widetilde{V}_0, V_{\alpha}]^{K+1}(x)s, \ldots, [\widetilde{V}_0, V_{\alpha}]^N(x)s ). \end{align*}$$ Under Assumption 1 it holds that for any |$\varepsilon \in [0,1]$|⁠, |$x \in \mathbf{R}^N$| and |$s>0$|⁠, $$\begin{equation*}c_{1,1}^{\varepsilon}(s,x), \ldots, c_{1,d}^{\varepsilon}(s,x), c_{2,1}^{\varepsilon}(s,x), \ldots, c_{2,d}^{\varepsilon}(s,x)\end{equation*}$$ linearly span |$\mathbf{R}^N$|⁠, which implies $$\begin{align*} \inf_{\|\xi\|=1, \xi \in \mathbf{R}^N} \sum_{j=1}^2 \sum_{\alpha =1}^d \langle c_{j, \alpha}^{\varepsilon}(s,x), \xi \rangle^2>0, \end{align*}$$ where |$\langle \cdot , \cdot \rangle $| is the Euclidean inner product. Hence, applying the result of Proposition 9 in Bally (2003) leads to the fact: for |$p<\infty $| there is |$C>0$| such that |$\sup _{\varepsilon \in [0,1]} \| (\det \sigma ^{Y_1^{\varepsilon }})^{-1} \|_p \leqslant C$| for all |$x \in \mathbf{R}^N$|⁠. Since |$\{ Y_1^{\epsilon }, \epsilon \in (0,1) \}$| is uniformly nondegenerate with respect to |$\epsilon $| we have the expansion of |$E[\delta _y (Y_1^{\varepsilon })]$| around |$E[\delta _y (Y_1^{0})]$| as follows: $$\begin{align} &E[\delta_y (Y_1^{\varepsilon})]=E[\delta_y (Y_1^{0})]+\sum_{j=1}^5 \varepsilon^{\,j} \sum_{k=1}^{\,j} \sum_{ \substack{ (\beta_i)_{i=1}^k \subset \mathbf{N} \\ s.t. \ \sum_{i=1}^k \beta_i=j}} \sum_{ (\alpha_1,\cdots,\alpha_k) \in \{1,\cdots,N \}^k } \frac{1}{k!} E\!\left[ \!\delta_y(Y_1^0) H_{(\alpha_1,\cdots,\alpha_k)}\!\left(\!Y_1^0,\prod_{l=1}^k\! Y_1^{\beta_l,\alpha_l}\!\right) \!\right] \nonumber \\ & \qquad\qquad\qquad+\varepsilon^{6} \int_0^1 (1-\eta)^{5} \sum_{ \substack{ k=1,\cdots,6, \\ \beta=(\beta_1,\cdots,\beta_k), \\ \alpha=(\alpha_1,\cdots,\alpha_k) }}^{(6)} E \left[ \delta_y(Y_1^{\varepsilon \eta}) H_{(\alpha_1,\cdots, \alpha_k)} \left( Y_1^{\varepsilon \eta},\prod_{l=1}^k G^{\varepsilon \eta,l} \right) \right] \textrm{d}\eta, \end{align}$$(A.5) where $$\begin{eqnarray*} \sum_{ \substack{ k=1,\cdots,m, \\ \beta=(\beta_1,\cdots,\beta_k), \\ \alpha=(\alpha_1,\cdots, \alpha_k)}}^{(m)}&\equiv& m \sum_{k=1}^{m} \sum_{ \substack{ (\beta_i)_{i=1}^k \subset \mathbf{N}, \\ s.t. \ \sum_{i=1}^k \beta_i=m}} \ \sum_{ (\alpha_1,\cdots,\alpha_k) \in \{1,\cdots,N \}^k } \frac{1}{k!}, m \in \mathbf{N}, \end{eqnarray*}$$|$Y_1^{\beta _l,\alpha _l}$| is given by (2.21) or (2.22) and |$G^{\varepsilon ,l}=\frac{1}{\beta _{l}!} \frac{\partial ^{\beta _l}}{\partial \varepsilon ^{\beta _l}} Y_1^{\alpha _l, \varepsilon }$|⁠, |$\varepsilon \in [0,1]$| satisfies that for |$k=1,\cdots ,6$|⁠, |$s \in \mathbf{N}$|⁠, |$p\geqslant 1$|⁠, there is |$C>0$| such that |$\big \| \prod _{l=1}^k G^{\varepsilon ,l} \big \|_{s,p}\leqslant C$|⁠. By (A.3) we have the expansion $$\begin{align*} &E[\varphi(X_1^{x,\varepsilon})]=E[\varphi(\bar{X}_1^{x,\varepsilon})]\!+\!\sum_{j=1}^5 \varepsilon^{\,j}\! \sum_{k=1}^{\,j} \sum_{ \substack{ (\beta_i)_{i=1}^k \subset \mathbf{N} \\ s.t. \ \sum_{i=1}^k \beta_i=j} } \sum_{ (\alpha_1,\cdots,\alpha_k) \in \{\!1,\cdots,N \}^k } \frac{1}{k!} E\!\left[ \!\varphi(\bar{X}_1^{x,\varepsilon}) H_{(\alpha_1,\cdots,\alpha_k)}\!\left(\!\!Y_1^{0},\!\prod_{l=1}^k Y_1^{\beta_l,\alpha_l}\!\!\right)\! \right]\\ & +\varepsilon^6 \int_0^1 (1-\eta)^{5} \sum_{ \substack{ k=1,\cdots,6, \\ \beta=(\beta_1,\cdots,\beta_k), \\ \alpha=(\alpha_1,\cdots,\alpha_k) }}^{(6)} E \left[ \varphi(\widetilde{X}_1^{x, \varepsilon,\eta}) H_{(\alpha_1,\cdots, \alpha_k)} \left( Y_1^{\varepsilon \eta},\prod_{l=1}^k G^{\varepsilon \eta,l} \right) \right] \textrm{d}\eta, \end{align*}$$ where the |$N$|-dimensional random vector |$\widetilde{X}_1^{x, \varepsilon ,\eta }=(\widetilde{X}_1^{x,1, \varepsilon ,\eta },\cdots ,\widetilde{X}_1^{x,N,\varepsilon ,\eta }),\ \epsilon , \eta \in [0,1]$| is given by |$\widetilde{X}_1^{x,i,\varepsilon ,\eta } = \tau ^i(Y_1^{\epsilon \eta }),\; i=1,\ldots , N$|⁠, with the function |$\tau $| defined by (2.23) and (2.24). For the error term the following holds: for |$p<\infty $| and |$\alpha \in \{1,\cdots ,N \}^k$|⁠, |$k\geqslant 1$| there is a |$C>0$| such that $$\begin{eqnarray*} \sup_{\lambda \in [0,1]} \left\| H_{(\alpha_1,\cdots, \alpha_k)} \left( Y_1^{\lambda},\prod_{l=1}^k G^{\lambda,l} \right) \right\|_p \leqslant C, \end{eqnarray*}$$ for all |$x \in \mathbf{R}^N$|⁠, where we used |$\| H_{(\alpha _1,\cdots , \alpha _k)} ( Y_1^{\lambda },\prod _{l=1}^k G^{\lambda ,l} ) \|_p \leqslant C \| (\det \sigma ^{Y_1^{\lambda }})^{-1} \|_s^{\gamma } $| for some |$C,s,\gamma>0$| for all |$\lambda \in [0,1]$|⁠, |$x \in \mathbf{R}^N$| with (A.4). Therefore, we obtain the assertion. |$\Box $| B. Proof of Lemma 2.5 First, we note $$\begin{align} I_{\beta} (1) = \int_{0 < t_{i_n +1} < \cdots < t_m < 1} I_{( \beta_{1}, \ldots, \beta_{i_n})} (t_{i_n +1} ) \textrm{d}t_{i_n +1} \cdots \textrm{d}t_{m} \nonumber \end{align}$$ due to the definition of |$\hat{\beta }$| given |$\beta \in\{0,1,\cdots ,d \}^m$|⁠. Then, by Fubini’s theorem, interchanging the order of integration with respect to the Lebesgue measure and the Wiener measure provides $$\begin{align} E [ \partial^{\gamma} \varphi ( \bar{X}_1^{x,\varepsilon} ) \varepsilon^{\| \beta \|} I_{\beta}(1) ] = \varepsilon^{\| \beta \|}\int_{0 < t_{i_n +1} < \cdots < t_m < 1} E [ \partial^{\gamma} \varphi ( \bar{X}_1^{x,\varepsilon} ) I_{( \beta_{1}, \ldots, \beta_{i_n})} (t_{i_n +1}) ] \textrm{d}t_{i_n +1} \cdots \textrm{d}t_{m}. \end{align}$$(B.1) Note that |$I_{( \beta _{1}, \ldots , \beta _{i_n})} (t_{i_n +1})$| is interpreted as a Skorohod integral, and hence the duality relation gives $$\begin{align} E [ \partial^{\gamma} \varphi( \bar{X}_1^{x,\varepsilon} ) I_{( \beta_{1}, \ldots, \beta_{i_n})} (t_{i_n +1}) ] = &\; E \left[ \int_0^{t_{i_n +1}} D_{t_{i_n}} ^{\beta_{i_n}} \partial^{\gamma} \varphi( \bar{X}_1^{x,\varepsilon} ) I_{( \beta_{1}, \ldots, \beta_{i_n-1})} (t_{i_n}) \textrm{d}t_{i_n} \right] \nonumber \\ =&\; \sum_{\alpha_n=1}^N E \left[ \partial_{\alpha_n} \partial^{\gamma} \varphi( \bar{X}_1^{x,\varepsilon} ) \int_0^{t_{i_n +1}} D_{t_{i_n}}^{\beta_{i_n}} \bar{X}_1^{x,\alpha_n,\varepsilon} (x) I_{( \beta_{1}, \ldots, \beta_{i_n-1})} (t_{i_n}) \textrm{d}t_{i_n} \right] \nonumber \\ =&\; \sum_{\alpha_n=1}^N E \left[ \partial_{\alpha_n} \partial^{\gamma} \varphi( \bar{X}_1^{x,\varepsilon} ) \int_0^{t_{i_n +1}} {\psi}_{\beta_{i_n}}^{\alpha_n, \varepsilon}(t_{i_n},x) I_{( \beta_{1}, \ldots, \beta_{i_n-1})} (t_{i_n}) \textrm{d}t_{i_n} \right], \end{align}$$(B.2) where |${ \psi }_{\beta _{i_n}}^{\alpha _n, \varepsilon }$| is the deterministic function given in Lemma 2.5. Repeating the same procedure for (B.2) until the multiple Itô integral is converted into the iterated time integral we obtain the formula (2.30). Note that when the multi-index |$\beta $| satisfies |$\# \{ j ; \beta _j=0 \}=i \leqslant m$|⁠, i.e. |$\| \beta \|=m+i$|⁠, |$m\geqslant 3$|⁠, the product |${{\psi }_{\beta _{i_1}}^{\alpha _1,\varepsilon }(t_{{i_1}},x) \cdots{\psi } _{\beta _{i_n}}^{\alpha _n,\varepsilon }(t_{i_n},x)}$| in (2.30) will be |$O(\varepsilon ^{m+\nu -i})$| for a |$\nu \geqslant 0$| by the definition of |${\psi }_{i_k}^{j,\varepsilon }(t,x)$|⁠. Then, when |$m\geqslant 3$| there exists a |$C>0$| such that $$\begin{eqnarray*} \sup_{x \in \mathbf{R}^N}\Big|E [ \partial^{\gamma} \varphi( \bar{X}_1^{x,\varepsilon} ) \varepsilon^{\| \beta \|} I_{\beta}(1) ] \Big| \leqslant C \varepsilon^{m+i+m+\nu-i} \leqslant C \varepsilon^{m+i+m-i}=C \varepsilon^{2m} \leqslant C \varepsilon^{6}, \end{eqnarray*}$$ for all |$\varepsilon \in (0,1]$|⁠, if we use the derivatives of |$\varphi $|⁠, and therefore we get the estimate (2.31). |$\Box $| C. Proof of Lemma 2.7 Let |$\{ \theta _t\}_{t \geqslant 0}$| be the stochastic process given by $$\begin{align} \theta_t = ( W_t^1, \ldots, W_t^d, {\int_0^t W_s^{\alpha_1} \textrm{d}s, \ldots, \ \int_0^t W_s^{\alpha_l} \textrm{d}s}), \end{align}$$(C.1) where |$\alpha _1, \ldots , \alpha _l$| are all elements in |$I = {\bigcup _{m=K+1}^N} \{ \alpha \in \{ 1, 2,\ldots , d \}: \ [V_{\alpha }, V_0 ]^m (\cdot ) \neq 0 \}$|⁠. Note that if |$\alpha \in I$|⁠, |${\int _0^t W_s^{\alpha } \textrm{d}s}$| appears in |$\bar{X}_t^x$|⁠, otherwise it does not. Let us denote the Malliavin covariance of |$\theta _t$| by |$\sigma ^{\theta _t} = ( \sigma ^{\theta _t}_{i\,j})_{1 \leqslant i,j \leqslant d+l}$|⁠. Then, one has its inverse denoted by |$\gamma ^{\theta _t} = ( \gamma ^{\theta _t}_{i\,j} )_{1 \leqslant i,j \leqslant d+l}$| as follows: $$\begin{align*} \gamma_{ii}^{\theta_t} = \left\{ \begin{array}{ll} {1}/{t} & \textrm{if} \ \ i \notin I \\{4}/{t} & \textrm{if} \ \ i \in I \\{12}/{t^3} & \textrm{if} \ \ i = d+1, \ldots, d+l \end{array} \right. \end{align*}$$ and |$\gamma ^{\theta _t}_{d + m\;i} = \gamma ^{\theta _t}_{i \;d+m} = -6/ t^2$| if |$i = \alpha _m,$||$m=1,2, \ldots , l$|⁠. For any other elements |$\gamma _{i\,j}^{\theta _t}$| takes |$0$|⁠. With the scaling property, we obtain $$\begin{align} &\; \sum_{i=1}^K E \left[ \partial_{i} \varphi ( \bar{X}_1^{x,\sqrt{t}} ) \sum_{0 \leqslant j_1,\,j_2 \leqslant d} t^{\frac{1}{2} \| (\,j_1,\,j_2)\|} \hat{V}_{j_1} V_{j_2}^{i}(x) I_{(j_1,j_2)}(1) \right] \nonumber\\ = &\; \sum_{i=1}^K E \left[ \partial_{i} \varphi ( \bar{X}_t^x ) \sum_{0 \leqslant j_1,\,j_2 \leqslant d} \hat{V}_{j_1} V_{j_2}^{i}(x) \int_0^t W_s^{j_1} \textrm{d}W_s^{j_2} \right] \nonumber \\ = &\; \sum_{i=1}^K \int_{\mathbf{R}^{d+l}} \partial_i \varphi \left( \zeta (y) \right) \; {\vphantom{\Bigg(}}_{-\infty} \left\langle \delta_y ( \theta_t ), \sum_{0 \leqslant j_1,\,j_2 \leqslant d} \hat{V}_{j_1} V_{j_2}^{i}(x) \int_0^t W_s^{j_1} \textrm{d}W_s^{j_2} \right\rangle_{ \infty} \textrm{d}y, \end{align}$$(C.2) where the function |$\zeta : \mathbf{R}^{d+l} \to \mathbf{R}^N$| is defined by $$\begin{align} \zeta_{i_1} (y) =&\; x_{i_1} + V_0^{i_1}(x)t + \sum_{j=1}^d V_j^{i_1} (x) y_{j}, \ \ i_1 =1,\ldots, K, \nonumber \\ \zeta_{i_2}(y) = &\; x_{i_2} + V_0^{i_2}(x)t + \hat{V}_0 V_0^{i_2}(x) \frac{t^2}{2} + \sum_{m=1}^l \hat{V}_{\alpha_m} V_0^{i_2} (x) y_{d+m}, \ \ i_2 = K+1, \ldots, N. \nonumber \end{align}$$ For the the coupling in (C.2) it holds $$\begin{align} &\; {\vphantom{\Bigg(}}_{-\infty} \left\langle \delta_y (\theta_t ), \sum_{0 \leqslant j_1,j_2 \leqslant d} \hat{V}_{j_1} V_{j_2}^{i}(x) \int_0^t W_s^{j_1} \textrm{d}W_s^{j_2} \right\rangle_{ \infty} \nonumber \\ &\quad= {\vphantom{\Bigg(}}_{-\infty} \left\langle \delta_y (\theta_t ), \sum_{0 \leqslant j_1,j_2 \leqslant d} \hat{V}_{j_1} V_{j_2}^{i}(x) \left\{ \frac{1}{2} W_t^{j_1} W_t^{j_2} - \frac{t}{2} \mathbf{1}_{j_1 = j_2 \neq 0} \right\} \right\rangle_{ \infty} \end{align}$$(C.3) $$\begin{align} &\; + \; {\vphantom{\Bigg(}}_{-\infty} \left\langle \delta_y (\theta_t ), \sum_{\alpha \in I} \sum_{\substack{ 0 \leqslant j \leqslant d \\ j \neq \alpha}} [V_{\alpha}, V_{j}]^i(x) \left\{ \frac{1}{t} I_{(\alpha,0)}(t) W_t^{j} - \frac{1}{2} W_t^{\alpha} W_t^{j} \right\} \right\rangle_{\infty}. \end{align}$$(C.4) Then, substituting this formula into (C.2), we easily obtain |$\pi _{1,t}^x$| in Lemma 2.4 by using integration by parts. Then, we will show that the coupling in the equation (C.2) is represented as a summation of (C.3) and (C.4). First of all note that we have, for |$G \in{{\mathbf D}}^{\infty }$|⁠, $$\begin{align} H_{(j)} (\theta_t, G) =&\; \frac{1}{t} G W_t^{\,j} - \frac{1}{t} \int_{0}^t D^{\,j}_s G \textrm{d}s, \ \ j \notin I, \end{align}$$(C.5) $$\begin{align} H_{(j)} (\theta_t, G) =&\; \frac{4}{t} G W_t^{\,j} - \frac{6}{t^2} G I_{(j,0)}(t) - \frac{4}{t} \int_{0}^t D^{\,j}_s G \textrm{d}s, + \frac{6}{t^2} \int_{0}^t D^{\,j}_s G (t-s) \textrm{d}s, \ \ j \in I, \end{align}$$(C.6) and for |$1 \leqslant m \leqslant l,$| $$\begin{align} H_{(d+m)} ( \theta_t, G ) = - \frac{6}{t^2}G W_t^{\alpha_m} + \frac{12}{t^3} G I_{(\alpha_m,0)} (t) + \frac{6}{t^2} \int_{0}^t D^{\alpha_m}_s G \textrm{d}s - \frac{12}{t^3} \int_0^t D_s^{\alpha_m} G (t-s) \textrm{d}s, \end{align}$$(C.7) where |$\alpha _m \in I$|⁠. Since we can derive the equations (C.5), (C.6) and (C.7) in almost the same way we will just give a computation for (C.7). One has $$\begin{eqnarray*} H_{(d+m)} ( \theta_t, G ) &=& \sum_{j=1}^{d+l} \delta ( \gamma_{d+m\; j}^{\theta_t} G D\theta_t^{\,j} ) = - \frac{6}{t^2} \delta ( G D\theta_t^{\alpha_m} ) + \frac{12}{t^3} \delta ( G D \theta_t^{d + m} ) \nonumber \\ &=& - \frac{6}{t^2} \left( G W_t^{\alpha_m} - \int_0^t D_s^{\alpha_m} G \textrm{d}s \right) + \frac{12}{t^3} \left( G I_{(\alpha_m,0)}(t) - \int_0^t D_s^{\alpha_m} G (t-s) \textrm{d}s \right), \nonumber \end{eqnarray*}$$ where we used the relation |$\delta (F h) = F W(h) - \langle DF, h \rangle _H$| for |$h \in H$| and |$F \in{{\mathbf D}}^{\infty }$|⁠. Now, the coupling in (C.2) are divided into four sums: $$\begin{equation} _{- \infty} \langle \delta_y (\theta_t ), \sum_{0 \leqslant j_1,\,j_2 \leqslant d} \hat{V}_{j_1} V_{j_2}^{i}(x) \int_0^t W_s^{j_1} \textrm{d}W_s^{j_2} \rangle_{ \infty} \end{equation}$$(C.8) $$\begin{align} =&\; \{ \sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1,\, j_2 \notin I}} + \sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1 \in I,\, j_2 \notin I}} + \sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1 \notin I,\, j_2 \in I}} + \sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1, \, j_2 \in I}} \} {\vphantom{\Bigg(}}_{-\infty} \left\langle \delta_y (\theta_t ), \hat{V}_{j_1} V_{j_2}^{i}(x) \int_0^t W_s^{j_1} \textrm{d}W_s^{j_2} \right\rangle_{ \infty}. \end{align}$$(C.9) For the first sum, we have $$\begin{align} &\sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1,\, j_2 \notin I}} {\vphantom{\Bigg(}}_{-\infty} \left\langle \delta_y (\theta_t ), \hat{V}_{j_1} V_{j_2}^{i}(x) \int_0^t W_s^{j_1} \textrm{d}W_s^{j_2} \right\rangle_{ \infty} = \sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1,\, j_2 \notin I}} {\vphantom{\Bigg(}}_{-\infty} \left\langle \partial_{j_2} \delta_y (\theta_t ), \hat{V}_{j_1} V_{j_2}^{i}(x) \int_0^t W_s^{j_1} \textrm{d}s \right\rangle_{ \infty} \nonumber \\ &= \sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1,\, j_2 \notin I}} {\vphantom{\Bigg(}}_{-\infty} \left\langle \partial_{j_2} \delta_y (\theta_t ), \hat{V}_{j_1} V_{j_2}^{i}(x) \int_0^t (t-s) \textrm{d}W_s^{j_1} \right\rangle_{ \infty} = \sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1,\, j_2 \notin I}} {\vphantom{\Bigg(}}_{-\infty} \left\langle \partial_{j_1} \partial_{j_2} \delta_y (\theta_t ), \hat{V}_{j_1} V_{j_2}^{i}(x) \frac{t^2}{2} \right\rangle_{ \infty}. \end{align}$$(C.10) In particular, we suppose for |$j =0$|⁠, |$\partial _j$| is an identity operator on |$\delta _y.$| In addition, since we have |$H_{(j_1,j_2)} ( \theta _t,1) = {\frac{1}{t^2}W_t^{j_1} W_t^{j_2} - \frac{1}{t} \mathbf{1}_{ j_1 = j_2 \neq 0}} $| for |$j_1,j_2 \notin I$|⁠, we get $$\begin{align} \sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1, \,j_2 \notin I}} {\vphantom{\Bigg(}}_{-\infty}\! \left\langle \!\delta_y (\theta_t ), \hat{V}_{j_1}\! V_{j_2}^{i}(x)\! \int_0^t \!W_s^{j_1} \textrm{d}W_s^{j_2} \right\rangle_{ \infty} \!=\! \sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1,\, j_2 \notin I}} {\vphantom{\Bigg(}}_{-\infty}\! \left\langle \delta_y (\theta_t ), \hat{V}_{j_1} \!V_{j_2}^{i}(x) \!\left\{ \frac{1}{2}W_t^{j_1} W_t^{j_2} \!-\! \frac{t}{2} \mathbf{1}_{ j_1 = j_2 \neq 0} \right\}\!\right\rangle_{ \infty}. \end{align}$$(C.11) For the second and third sum of (C.9) we compute as follows: $$\begin{align} &\; \left\{ \sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1 \in I,\, j_2 \notin I}} + \sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1 \notin I, j_2 \in I}} \right\} \; {\vphantom{\Bigg(}}_{-\infty} \left\langle \delta_y (\theta_t ), \hat{V}_{j_1} V_{j_2}^{i}(x) \int_0^t W_s^{j_1} \textrm{d}W_s^{j_2} \right\rangle_{ \infty} \nonumber \\ =&\; \sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1 \in I,\, j_2 \notin I}} \; {\vphantom{\Bigg(}}_{-\infty} \left\langle \delta_y (\theta_t ), \hat{V}_{j_1} V_{j_2}^{i}(x) \int_0^t W_s^{j_1} \textrm{d}W_s^{j_2} \right\rangle_{ \infty} \nonumber\\&+ \sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1 \notin I,\, j_2 \in I}} \; {\vphantom{\Bigg(}}_{-\infty} \left\langle \delta_y (\theta_t ), \hat{V}_{j_1} V_{j_2}^{i}(x) \left\{ W_t^{j_1} W_t^{j_2} - \int_0^t W_s^{j_2} \textrm{d}W_s^{j_1} \right\} \right\rangle_{ \infty} \nonumber \\ =&\; \sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1 \in I,\, j_2 \notin I}} \; {\vphantom{\Bigg(}}_{-\infty} \left\langle \delta_y (\theta_t ), [V_{j_1}, V_{j_2} ]^i (x) \int_0^t W_s^{j_1} \textrm{d}W_s^{j_2} \right\rangle_{ \infty} \nonumber\\&+ \sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1 \in I,\, j_2 \notin I}} {}_{- \infty} \langle \delta_y (\theta_t ), \hat{V}_{j_2} V_{j_1}^i(x) W_t^{j_1} W_t^{j_2} \rangle_{ \infty}. \end{align}$$(C.12) Applying the integration by parts to the first term of (C.12), we get $$\begin{align} \sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1 \in I,\, j_2 \notin I}} &\; {\vphantom{\Bigg(}}_{-\infty} \left\langle \delta_y (\theta_t ), [V_{j_1}, V_{j_2} ]^i (x) \int_0^t W_s^{j_1} \textrm{d}W_s^{j_2} \right\rangle_{ \infty}\nonumber\\ &= \sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1 \in I,\, j_2 \notin I}} \; {\vphantom{\Bigg(}}_{-\infty} \left\langle \partial_{j_2} \delta_y (\theta_t ), [V_{j_1}, V_{j_2} ]^i (x) \int_0^t W_s^{j_1} \textrm{d}s \right\rangle_{ \infty} \nonumber \\ &=\sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1 \in I,\, j_2 \notin I}} \; {\vphantom{\Bigg(}}_{-\infty} \left\langle \delta_y (\theta_t ), [V_{j_1}, V_{j_2} ]^i (x) H_{(j_2)} \left(\theta_t, I_{(j_1,0)}(t) \right) \right\rangle_{ \infty}.\end{align}$$(C.13) Since we have |$ H_{(j_2)} \left ( \theta _t, I_{(j_1,0)}(t) \right ) = { \frac{1}{t} I_{(j_1,0)}(t) W_t^{j_2}} $| for |$j_1 \in I$| and |$j_2 \notin I$| we deduce that $$\begin{align} & \left\{ \sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1 \in I, \,j_2 \notin I}} + \sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1 \notin I, j_2 \in I}} \right\} \; {\vphantom{\Bigg(}}_{-\infty} \left\langle \delta_y (\theta_t ), \hat{V}_{j_1} V_{j_2}^{i}(x) \int_0^t W_s^{j_1} \textrm{d}W_s^{j_2} \right\rangle_{ \infty} \nonumber \\ =&\; \sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1 \in I,\, j_2 \notin I}} \; {\vphantom{\Bigg(}}_{-\infty} \left\langle \delta_y (\theta_t ), [V_{j_1}, V_{j_2} ]^i (x) \frac{1}{t} I_{(j_1,0)}(t) W_t^{j_2} \right\rangle_{ \infty} + \sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1 \in I,\, j_2 \notin I}} {}_{- \infty} \langle \delta_y (\theta_t ), \hat{V}_{j_2} V_{j_1}^i(x) W_t^{j_1} W_t^{j_2} \rangle_{ \infty} \nonumber \\ =&\; \sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1 \in I,\, j_2 \notin I}} \; {\vphantom{\Bigg(}}_{-\infty} \left\langle \delta_y (\theta_t ), [V_{j_1}, V_{j_2} ]^i (x) \left\{ \frac{1}{t} I_{(j_1,0)}(t) W_t^{j_2} - \frac{1}{2}W_t^{j_1} W_t^{j_2} \right\} \right\rangle_{ \infty} \nonumber \\ &\; + \sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1 \in I,\, j_2 \notin I}} \; {\vphantom{\Bigg(}}_{-\infty} \left\langle \delta_y (\theta_t ), [V_{j_1}, V_{j_2} ]^i (x) \frac{1}{2}W_t^{j_1} W_t^{j_2} \right\rangle_{ \infty} + \sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1 \in I,\, j_2 \notin I}} {}_{- \infty} \langle \delta_y (\theta_t ), \hat{V}_{j_2} V_{j_1}^i(x) W_t^{j_1} W_t^{j_2} \rangle_{ \infty} \nonumber \\ =&\; \sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1 \in I, \,j_2 \notin I}} \; {\vphantom{\Bigg(}}_{-\infty} \left\langle \delta_y (\theta_t ), [V_{j_1}, V_{j_2} ]^i (x) \left\{ \frac{1}{t} I_{(j_1,0)}(t) W_t^{j_2} - \frac{1}{2}W_t^{j_1} W_t^{j_2} \right\} \right\rangle_{ \infty} \nonumber\\ &\; +\left\{ \sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1 \in I,\, j_2 \notin I}} + \sum_{\substack{0 \leqslant j_1,\,j_2 \leqslant d \\ j_1 \notin I,\, j_2 \in I}} \right\} \; {\vphantom{\Bigg(}}_{-\infty} \left\langle \delta_y (\theta_t ), \hat{V}_{j_1} V_{j_2}^i(x)\frac{1}{2} W_t^{j_1} W_t^{j_2} \right\rangle_{ \infty}. \end{align}$$(C.14) Finally, we move on the fourth sum of (C.9), which is divided into two parts as follows: $$\begin{align} & \sum_{\substack{ j_1, j_2 \in I}} \; {\vphantom{\Bigg(}}_{-\infty} \left\langle \delta_y (\theta_t ), \hat{V}_{j_1} V_{j_2}^{i}(x) \int_0^t W_s^{j_1} \textrm{d}W_s^{j_2} \right\rangle_{ \infty} \nonumber \\ =&\; \sum_{\substack{ j \in I}} \; {\vphantom{\Bigg(}}_{-\infty} \left\langle \delta_y (\theta_t ), \hat{V}_{j} V_{j}^{i}(x) \left\{ \frac{1}{2} (W_t^{\,j})^2 - \frac{t}{2} \right\} \right\rangle_{ \infty} + \sum_{\substack{ j_1 \neq j_2 \in I}} {\vphantom{\Bigg(}}_{-\infty} \left\langle \delta_y (\theta_t ), \hat{V}_{j_1} V_{j_2}^{i}(x) \int_0^t W_s^{j_1} \textrm{d}W_s^{j_2} \right\rangle_{ \infty}.\nonumber \end{align}$$ The first term is derived in the case where |$j_1 =j_2 \in I$|⁠. Hereafter, we denote |$j_1 \neq j_2 \in I$| by |$\alpha _m,\alpha _n \in I$| where |$ 1 \leqslant m \neq n \leqslant l$|⁠. Then, the second term is calculated in the following way. $$\begin{align} &\sum_{\substack{ 1 \leqslant m,n \leqslant l \\ m \neq n}} {\vphantom{\Bigg(}}_{-\infty} \left\langle \delta_y (\theta_t ), \hat{V}_{\alpha_m} V_{\alpha_n}^{i}(x) \int_0^t W_s^{\alpha_m} \textrm{d}W_s^{\alpha_n} \right\rangle_{ \infty} \nonumber \\ =&\; \sum_{\substack{1 \leqslant m, n \leqslant l \\ m \neq n}} \; {\vphantom{\Bigg(}}_{-\infty} \left\langle \partial_{\alpha_n} \delta_y (\theta_t ), \hat{V}_{\alpha_m} V_{\alpha_n}^{i}(x) \int_0^t W_s^{\alpha_m} \textrm{d}s \right\rangle_{ \infty} \nonumber\\&+ \sum_{\substack{1 \leqslant m, n \leqslant l \\ m \neq n}} \; {\vphantom{\Bigg(}}_{-\infty} \left\langle \partial_{d+n} \delta_y (\theta_t ), \hat{V}_{\alpha_m} V_{\alpha_n}^{i}(x) \int_0^t (t-s) W_s^{\alpha_m} \textrm{d}s \right\rangle_{ \infty} . \end{align}$$(C.15) Here, Itô formula gives |$\int _0^t (t-s) W_s^{\alpha _m} \textrm{d}s = \int _0^t \frac{1}{2}(t-s)^2 \textrm{d}W_s^{\alpha _m} $| and then use of integration by parts leads to $$\begin{align} &(\text{2nd term of (C.15)}) \ = \sum_{\substack{1 \leqslant m, n \leqslant l \\ m \neq n}} \; {\vphantom{\Bigg(}}_{-\infty} \left\langle \partial_{d+n} \delta_y (\theta_t ), \hat{V}_{\alpha_m} V_{\alpha_n}^{i}(x) \int_0^t \frac{1}{2} (t-s)^2 \textrm{d}W_s^{\alpha_m} \right\rangle_{ \infty} \nonumber \\ &=\; \sum_{\substack{1 \leqslant m, n \leqslant l \\ m \neq n}} \; {\vphantom{\Bigg(}}_{-\infty} \left\langle \partial_{\alpha_m} \partial_{d+n} \delta_y (\theta_t ), \hat{V}_{\alpha_m} V_{\alpha_n}^{i}(x) \frac{t^3}{6} \right\rangle_{ \infty} + \sum_{\substack{1 \leqslant m, n \leqslant l \\ m \neq n }} \; {\vphantom{\Bigg(}}_{-\infty} \left\langle \partial_{d+m} \partial_{d+n} \delta_y (\theta_t ), \hat{V}_{\alpha_m} V_{\alpha_n}^{i}(x) \frac{t^4}{8} \right\rangle_{ \infty}. \nonumber \end{align}$$ Hence, we obtain $$\begin{align} & \sum_{\substack{ 1 \leqslant m,n \leqslant l \\ m \neq n}} \; {\vphantom{\Bigg(}}_{-\infty}\left\langle \delta_y (\theta_t ), \hat{V}_{\alpha_m} V_{\alpha_n}^{i}(x) \int_0^t W_s^{\alpha_m} \textrm{d}W_s^{\alpha_n} \right\rangle_{ \infty} \nonumber \\ =&\; \sum_{\substack{ 1 \leqslant m,n \leqslant l \\ m \neq n}} \; {\vphantom{\Bigg(}}_{-\infty} \left\langle \delta_y (\theta_t ), \; \hat{V}_{\alpha_m} V_{\alpha_n}^{i}(x) \left\{ H_{( \alpha_n)}\left( \theta_t, I_{(\alpha_m,0)}(t) \right) + H_{( \alpha_m, d+n)}( \theta_t, 1 )\frac{t^3}{6} + H_{( d+m, d+n)}( \theta_t, 1 )\frac{t^4}{8} \right\} \right\rangle_{ \infty} \nonumber \\ =&\; \sum_{\substack{ 1 \leqslant m,n \leqslant l \\ m \neq n}} \; {\vphantom{\Bigg(}}_{-\infty} \left\langle \delta_y (\theta_t ), \; \hat{V}_{\alpha_m} V_{\alpha_n}^{i}(x) \left\{ \frac{1}{2} W_t^{\alpha_m} W_t^{\alpha_n} + \frac{1}{t} I_{(\alpha_m,0)}(t) W_t^{\alpha_n} - \frac{1}{t} W_t^{\alpha_m} I_{(\alpha_n,0)}(t) \right\} \right\rangle_{ \infty} \nonumber \\ =&\; \sum_{\substack{ 1 \leqslant m,n \leqslant l \\ m \neq n}} \; {\vphantom{\Bigg(}}_{-\infty} \left\langle \delta_y (\theta_t ), \; \hat{V}_{\alpha_m} V_{\alpha_n}^{i}(x) \frac{1}{2} W_t^{\alpha_m} W_t^{\alpha_n} \right\rangle_{ \infty} \nonumber \\ &\;+ \sum_{\substack{ 1 \leqslant m,n \leqslant l \\ m \neq n}} \; {\vphantom{\Bigg(}}_{-\infty} \left\langle \delta_y (\theta_t ), [V_{\alpha_m}, V_{\alpha_n} ]^i (x) \left\{ \frac{1}{t} I_{(\alpha_m,0)}(t) W_t^{\alpha_n} - \frac{1}{2} W_t^{\alpha_m} W_t^{\alpha_n} \right\} \right\rangle_{ \infty}. \end{align}$$(C.16) On the second equation above we used the following results: $$\begin{eqnarray*} &&H_{(\alpha_n)} \left(\theta_t, I_{(\alpha_m,0)}(t) \right) = \frac{4}{t} I_{(\alpha_m,0)}(t) W_t^{\alpha_n} -\frac{6}{t^2}I_{(\alpha_m,0)}(t) I_{(\alpha_n,0)}(t), \nonumber \\ &&H_{(\alpha_m, d+n)} ( \theta_t,1 ) = - \frac{24}{t^3} W_t^{\alpha_m} W_t^{\alpha_n} + \frac{36}{t^4} I_{(\alpha_m,0)}(t)W_t^{\alpha_n} + \frac{48}{t^4} W_t^{\alpha_m}I_{(\alpha_n,0)}(t) -\frac{72}{t^5}I_{(\alpha_m,0)}(t)I_{(\alpha_n,0)}(t), \nonumber \\ &&H_{(d+m, d+n)} ( \theta_t,1 ) = \frac{36}{t^4} W_t^{\alpha_m} W_t^{\alpha_n} - \frac{72}{t^5} I_{(\alpha_m,0)}(t) W_t^{\alpha_n} - \frac{72}{t^5} W_t^{\alpha_m} I_{(\alpha_n,0)}(t) + \frac{144}{t^6}I_{(\alpha_m,0)}(t)I_{(\alpha_n,0)}(t). \nonumber \end{eqnarray*}$$ Combining the results (C.11), (C.14), (C.15) and (C.17) we finally get the appropriate transformation of the coupling in (C.2) as a summation of (C.3) and (C.4). We immediately obtain the Malliavin weight |$\pi _{2,t}^x$| from the following computation: $$\begin{align*} E [\partial_{i_1} \partial_{i_2} \varphi (\bar{X}_1^{x, \sqrt{t}}) ] = E [\partial_{i_1} \partial_{i_2} \varphi (\bar{X}_t^{x}) ] =E[ \varphi( \bar{X}_t^x) H_{(i_1, i_2)}( \bar{X}_t^x,1)]. \end{align*}$$ D. Proof of Lemma 2.9 We will show that there exist a constant |$\varepsilon>0$| and an integer |$\mu $| such that |${\inf _{\|z\|=1, z \in \mathbf{R}^N} \langle \sigma ^{\bar{X}_{kT/n}^{(n),x}}z, z\rangle }> \varepsilon n^{-\mu }$| a.s. for |$k = [n/2], \ldots , n-1$|⁠. At first, from (2.40) and (2.41), the Malliavin derivative of |$\bar{X}_{kT/n}^{(n),x}$| is given by $$\begin{align} D_s^l \bar{X}_{kT/n}^{ (n),x, i_1} & = V_l^{i_1}( \bar{X}_{(k-1)T/n}^{(n),x}) \mathbf{1}_{[(k-1)T/n,kT/n]}(s) + \varPsi_l^{i_1} \mathbf{1}_{[0, (k-1)T/n]}(s), \ \ i_1 =1, \ldots, K, \nonumber \\ D_s^l \bar{X}_{kT/n}^{(n), x, i_2} &= \hat{V}_l V_0^{i_2} ( \bar{X}_{(k-1)T/n}^{(n),x})(kT/n -s) \mathbf{1}_{[(k-1)T/n,kT/n]}(s) + \varPsi_l^{i_2} \mathbf{1}_{[0, (k-1)T/n]}(s), \ \ i_2 =K+1, \ldots, N, \nonumber \end{align}$$ where |$\varPsi _l^{i}$| is explicitly given as the Wiener functional with first-order partial derivatives of |$V_0, \ldots , V_d$|⁠, |$\hat{V}_0 V_0, \ldots , \hat{V}_d V_0 $| depending on |$\bar{X}_{(k-1)T/n}^{(n),x}$| and its Malliavin derivative. Then, the Malliavin covariance of |$\bar{X}_{kT/n}^{(n),x}$| is of the form |$\sigma ^{\bar{X}_{kT/n}^{(n),x}} = A + B $|⁠, where the symmetric matrices |$A=(A_{i\,j})_{1 \leqslant i,j \leqslant N}$| and |$B =(B_{i\,j})_{1 \leqslant i,j \leqslant N}$| are given by $$\begin{equation} A_{i\,j} = \left\{ \begin{array}{lll} \sum_{l=1}^d V_l^i ( \bar{X}_{(k-1)T/n}^{(n),x}) V_l^{\,j} ( \bar{X}_{(k-1)T/n}^{(n),x}) \frac{T}{n} & 1 \leqslant i,j \leqslant K \\ \sum_{l=1}^d V_l^i ( \bar{X}_{(k-1)T/n}^{(n),x}) \hat{V}_l V_0^{\,j} ( \bar{X}_{(k-1)T/n}^{(n),x})\frac{1}{2} \left( \frac{T}{n} \right)^2 & 1 \leqslant i \leqslant K, \; K+1 \leqslant j \leqslant N \\ \sum_{l=1}^d \hat{V}_l V_0^i ( \bar{X}_{(k-1)T/n}^{(n),x}) \hat{V}_l V_0^{\,j} ( \bar{X}_{(k-1)T/n}^{(n),x}) \frac{1}{3} \left( \frac{T}{n} \right)^3 & K+1 \leqslant i,j \leqslant N, \end{array} \right. \end{equation}$$(D.1) and |$B_{i\,j} = \langle \varPsi ^{i} \mathbf{1}_{[0, (k-1)T/n]}, \varPsi ^{j} \mathbf{1}_{[0, (k-1)T/n]} \rangle _H$| for |$1 \leqslant i,j \leqslant N$|⁠. Let us denote by |$I_m$| the identity matrix of size |$m \times m$| and for any matrices |$A_1$| and |$A_2$| with the same size, we write |$A_1 \geqslant A_2$| if and only if |$A_1 - A_2$| is non-negative definite. Then, Assumption 1 guarantees that we can take some |$\varepsilon _1>0$| and |$\varepsilon _2>0$| for all |$x \in \mathbf{R}^N$| such that |$\sum _{k=1}^d V_k (x) \otimes V_k(x) \geqslant \varepsilon _1 I_{K}$| and |$\sum _{k=1}^d \bar{V}_k(x) \otimes \bar{V}_k(x) \geqslant \varepsilon _2 I_{N-K}$|⁠, where |$\bar{V}_j (x) = ( \hat{V}_j V_0^{K+1}(x), \ldots , \hat{V}_j V_0^N(x) ),\; j=1,\ldots ,d$|⁠. Hence, there exists |$\varepsilon>0$| such that $$\begin{align} A \geqslant \varepsilon \frac{1}{n^3} I_N \ \ \text{a.s.} \end{align}$$(D.2) Since the matrix |$B$| is symmetric non-negative definite by its definition, we obtain $$\begin{align} \inf_{ \| z \| =1, z \in \mathbf{R}^N} \langle \sigma^{\bar{X}_{kT/n}^{(n),x}} z, z \rangle \geqslant \inf_{ \| z \| =1, z \in \mathbf{R}^N} \langle A z, z \rangle \geqslant \varepsilon \frac{1}{n^3} \ \ \text{a.s.} \end{align}$$(D.3) This implies that the smallest eigenvalue of |$\sigma ^{\bar{X}_{kT/n}^{(n),x}}$| is no less than |$\varepsilon /n^3$| and then we easily deduce (2.46). The nondegeneracy of |${X}_{t}^{x}$|⁠, |$T/2 \leqslant t \leqslant T$| and the bound (2.47) follow from Kusuoka & Stroock (1985). E. Proof of Proposition 2.12 We note that for |$h \in H$| and |$F \in \mathbf{D}^{1,2}$| it follows that $$\begin{align} \delta (F h) = F W(h) - \langle DF, h \rangle_H. \end{align}$$(E.1) Then, using the formula (C.1), we obtain $$\begin{align*} & H_{(i)} \Big( \bar{X}_t^x, \frac{1}{2}{W_t^{j_1}W_t^{j_2}}- \frac{t}{2} \mathbf{1}_{ j_1 = j_2 \neq 0} \Big) \\ =& \sum_{l=1}^N \delta \left( \gamma^{\bar{X}_t^x}_{il} \Big( \frac{1}{2}{W_t^{j_1}W_t^{j_2}}- \frac{t}{2} \mathbf{1}_{ j_1 = j_2 \neq 0} \Big) D \bar{X}_t^{x, l} \right) \\ =& \sum_{l=1}^N \gamma^{\bar{X}_t^x}_{il} \left( \Big( \frac{1}{2}{W_t^{j_1}W_t^{j_2}}- \frac{t}{2} \mathbf{1}_{ j_1 = j_2 \neq 0} \Big) \int_0^t D_s \bar{X}_t^{x,l} \cdot \textrm{d}W_s - \frac{ W_t^{j_1}} {2} \int _0^t D_s^{j_2} \bar{X}_t^{x,l} \textrm{d}s - \frac{ W_t^{j_2}} {2} \int _0^t D_s^{j_1} \bar{X}_t^{x, l} \textrm{d}s \right) \\ =&\; \sum_{l=1}^K \gamma^{\bar{X}_t^x}_{il} \left( \Big( \frac{1}{2}{W_t^{j_1}W_t^{j_2}}- \frac{t}{2} \mathbf{1}_{ j_1 = j_2 \neq 0} \Big) \sum_{m=1}^d V_m^l (x) W_t^m - \frac{1}{2} V_{j_2}^l(x) W_t^{j_1} t - \frac{1}{2} V_{j_1}^l(x) W_t^{j_2} t \right) \\ & + \sum_{l=K + 1}^N \gamma^{\bar{X}_t^x}_{il} \left(\! \Big( \frac{1}{2}{W_t^{j_1}W_t^{j_2}}- \frac{t}{2} \mathbf{1}_{ j_1 = j_2 \neq 0} \Big)\! \sum_{m=1}^d \hat{V}_m V_0^l(x)\! \int_0^t W_s^m \textrm{d}s - \frac{1}{4} \hat{V}_{j_2} V_0^l (x) W_t^{j_1} t^2 \!-\! \frac{1}{4} \hat{V}_{j_1} V_0^l (x) W_t^{j_2} t^2 \right) \\ =&\; \Big( \frac{1}{2}{W_t^{j_1}W_t^{j_2}}- \frac{t}{2} \mathbf{1}_{ j_1 = j_2 \neq 0} \Big) \sum_{m=1}^d \left\{ \sum_{l=1}^K \gamma^{\bar{X}_t^x}_{il} V_m^l (x) W_t^m + \sum_{l=K + 1}^N \gamma^{\bar{X}_t^x}_{il} \hat{V}_m V_0^l(x) \int_0^t W_s^m \textrm{d}s \right\} \\ & - \frac{1}{2} \sum_{l=1}^K \gamma^{\bar{X}_t^x}_{il} \left( V_{j_2}^l(x) W_t^{j_1} + V_{j_1}^l(x) W_t^{j_2} \right) t -\frac{1}{4} \sum_{l=K + 1}^N \gamma^{\bar{X}_t^x}_{il} \left( \hat{V}_{j_2} V_0^l (x) W_t^{j_1} + \hat{V}_{j_1} V_0^l (x) W_t^{j_2} \right) t^2. \end{align*}$$ Next, we show the formula (2.54). First, we have $$\begin{align*} \;H_{(i_1)} ( \bar{X}_t^x,1 ) = &\; \sum_{l=1}^N \delta ( \gamma^{\bar{X}_t^x}_{i_1l} D \bar{X}_t^{x, l} ) \\ =&\; \sum_{l=1}^K \sum_{m=1}^d \gamma^{\bar{X}_t^x}_{i_1l}V_m^l (x) W_t^m + \sum_{l=K+1}^N \sum_{m=1}^d \gamma^{\bar{X}_t^x}_{i_1l} \hat{V}_m V_0^l(x) \int_0^t W_s^m \textrm{d}s. \end{align*}$$ Then, we get $$\begin{align*} H_{(i_1,i_2)}( \bar{X}_t^x, 1) = &\; H_{(i_2)} ( \bar{X}_t^x, H_{(i_1)} (\bar{X}_t^x,1 )) \\ = & \; \sum_{l=1}^N \delta \left( \gamma^{\bar{X}_t^x}_{i_2l} H_{(i_1)} (\bar{X}_t^x,1 ) D \bar{X}_t^{x, l} \right) \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \end{align*}$$ $$\begin{align*} \quad\quad\quad\quad\quad\quad\quad\quad=&\; H_{(i_1)} (\bar{X}_t^x,1 ) \sum_{l=1}^N \delta \left( \gamma^{\bar{X}_t^x}_{i_2l} D \bar{X}_t^{x, l} \right) -\sum_{l_1=1}^N \gamma^{\bar{X}_t^x}_{i_2 l_1} \sum_{m=1}^d \int_0^t D_s^{m} H_{(i_1)} (\bar{X}_t^x,1 ) \; D_s^{m} \bar{X}_t^{x, l_1} \textrm{d}s \\ =&\; H_{(i_1)} (\bar{X}_t^x,1 ) H_{(i_2)} (\bar{X}_t^x,1 ) -\sum_{l_1, l_2 =1}^N \gamma^{\bar{X}_t^x}_{i_2l_1} \gamma^{\bar{X}_t^x}_{i_1l_2} \sum_{m=1}^d \int_0^t D_s^m \bar{X}_t^{x, l_2} D_s^{m} \bar{X}_t^{x, l_1} \textrm{d}s \\ =&\; H_{(i_1)} (\bar{X}_t^x,1 ) H_{(i_2)} (\bar{X}_t^x,1 ) -\sum_{l_1, l_2 =1}^N \gamma^{\bar{X}_t^x}_{i_2l_1} \gamma^{\bar{X}_t^x}_{i_1l_2} \sigma_{l_1 l_2}^{\bar{X}_t^x} \\ =&\; H_{(i_1)} (\bar{X}_t^x,1 ) H_{(i_2)} (\bar{X}_t^x,1 ) - \gamma^{\bar{X}_t^x}_{i_1 i_2}. \end{align*}$$ Finally, we prove the formula (2.55). Similarly to the above, we have $$\begin{align*} & H_{(i)} \left( \bar{X}_t^x, \frac{1}{t} \left(\int_0^t W_s^{\alpha} \textrm{d}s \right) W_t^{j}- \frac{1}{2} W_t^{\alpha} W_t^{\,j} \right) \\ =&\; \sum_{l=1}^N \delta \left( \gamma^{\bar{X}_t^x}_{i l} \Big\{ \frac{1}{t} \Big(\int_0^t W_s^{\alpha} \textrm{d}s \Big) W_t^{j}- \frac{1}{2} W_t^{\alpha} W_t^{\,j} \Big\} D \bar{X}_t^{x, l} \right) \\ = &\; \sum_{l=1}^N \gamma^{\bar{X}_t^x}_{i l} \Big( \frac{1}{t} \Big(\int_0^t W_s^{\alpha} \textrm{d}s \Big) W_t^{j}- \frac{1}{2} W_t^{\alpha} W_t^{\,j} \Big) \int_0^t D_s \bar{X}_t^{x, l} \cdot \textrm{d}W_s \\& - \sum_{l=1}^N \gamma^{\bar{X}_t^x}_{i l} \left\{ \Big( \frac{1}{t} \int_0^t (t-s) D_s^{\alpha} \bar{X}_t^{x, l} \textrm{d}s -\frac{1}{2} \int_0^t D_s^{\alpha} \bar{X}_t^{x, l} \textrm{d}s \Big) W_t^{j} + \Big( \frac{1}{t} \Big(\int_0^t W_s^{\alpha} \textrm{d}s \Big) - \frac{1}{2} W_t^{\alpha} \Big) \int_0^t D_s^{\,j} \bar{X}_t^{x, l} \textrm{d}s \right\} \\ = &\; \Big( \frac{1}{t} \Big(\int_0^t W_s^{\alpha} \textrm{d}s \Big) W_t^{j}- \frac{1}{2} W_t^{\alpha} W_t^{\,j} \Big) \sum_{m=1}^d \left\{ \sum_{l=1}^K \gamma^{\bar{X}_t^x}_{i l} V_m^l (x) W_t^m + \sum_{l=K+1}^N \gamma^{\bar{X}_t^x}_{i l} \hat{V}_m V_0^l (x) \int_0^t W_s^m \textrm{d}s \right\} \\ & - \left( \frac{1}{t} \Big(\int_0^t W_s^{\alpha} ds \Big) - \frac{1}{2} W_t^{\alpha} \right) \sum_{l=1}^K \gamma^{\bar{X}_t^x}_{i l} V_j^l (x) t - \sum_{l=K+1}^N \gamma^{\bar{X}_t^x}_{i l} \frac{1}{12} \hat{V}_\alpha V_0^l (x) W_t^{\,j} t^2 \\ & - \left( \frac{1}{t} \Big(\int_0^t W_s^{\alpha} \textrm{d}s \Big) - \frac{1}{2} W_t^{\alpha} \right) \sum_{l= K+1}^N \gamma^{\bar{X}_t^x}_{i l} \frac{1}{2} \hat{V}_j V_0^l(x) t^2 \\ = &\; \Big( \frac{1}{t} \Big(\int_0^t W_s^{\alpha} \textrm{d}s \Big) W_t^{j}- \frac{1}{2} W_t^{\alpha} W_t^{\,j} \Big) \sum_{m=1}^d \left\{ \sum_{l=1}^K \gamma^{\bar{X}_t^x}_{i l} V_m^l (x) W_t^m + \sum_{l=K+1}^N \gamma^{\bar{X}_t^x}_{i l} \hat{V}_m V_0^l (x) \int_0^t W_s^m \textrm{d}s \right\} \\ & - \left(\! \frac{1}{t} \Big(\!\int_0^t W_s^{\alpha} \textrm{d}s \!\Big) \!-\! \frac{1}{2} W_t^{\alpha} \!\right)\! \left\{ \sum_{l=1}^K \gamma^{\bar{X}_t^x}_{i l} V_j^l (x) t \!+\! \sum_{l= K+1}^N \!\gamma^{\bar{X}_t^x}_{i l} \frac{1}{2} \hat{V}_j V_0^l(x) t^2 \right\} \!-\! \sum_{l=K+1}^N \!\gamma^{\bar{X}_t^x}_{i l} \frac{1}{12} \hat{V}_\alpha V_0^l (x) W_t^{\,j} t^2. \ \ \ \end{align*}$$ © The Author(s) 2020. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - A second-order discretization for degenerate systems of stochastic differential equations
JF - IMA Journal of Numerical Analysis
DO - 10.1093/imanum/draa039
DA - 2020-08-06
UR - https://www.deepdyve.com/lp/oxford-university-press/a-second-order-discretization-for-degenerate-systems-of-stochastic-bl9E98oidC
SP - 1
EP - 1
VL - Advance Article
IS - 
DP - DeepDyve
ER -