remarkRemark hypothesisHypothesis claimClaim High order splitting methods for
commutative SDEsJ. Foster and G. dos Reis and C. Strange

# High order splitting methods for SDEs
satisfying a commutativity condition††thanks:  Funding: The first author was
supported by the Department of Mathematical Sciences at the University of Bath
as well as the DataSig programme under the ESPRC grant EP/S026347/1. The
second author acknowledges support from the _Fundaç $\tilde{\text{a}}$o para a
Ci$\hat{e}$ncia e a Tecnologia_ (Portuguese Foundation for Science and
Technology) through the projects UIDB/00297/2020 and UIDP/00297/2020 (Center
for Mathematics and Applications, CMA/FCT/UNL).

James Foster .  University of Bath, Department of Mathematical Sciences.
jmf68@ bath.ac.uk Gonçalo dos Reis  University of Edinburgh, School of
Mathematics. G.dosReis@ ed.ac.uk Calum Strange  University of Edinburgh,
School of Mathematics. c.strange-1@ ed.ac.uk

###### Abstract

In this paper, we introduce a new simple approach to developing and
establishing the convergence of splitting methods for a large class of
stochastic differential equations (SDEs), including additive, diagonal and
scalar noise types. The central idea is to view the splitting method as a
replacement of the driving signal of an SDE, namely Brownian motion and time,
with a piecewise linear path that yields a sequence of ODEs – which can be
discretised to produce a numerical scheme. This new way of understanding
splitting methods is inspired by, but does not use, rough path theory. We show
that when the driving piecewise linear path matches certain iterated
stochastic integrals of Brownian motion, then a high order splitting method
can be obtained. We propose a general proof methodology for establishing the
strong convergence of these approximations that is akin to the general
framework of Milstein and Tretyakov. That is, once local error estimates are
obtained for the splitting method, then a global rate of convergence follows.
This approach can then be readily applied in future research on SDE splitting
methods. By incorporating recently developed approximations for iterated
integrals of Brownian motion into these piecewise linear paths, we propose
several high order splitting methods for SDEs satisfying a certain
commutativity condition. In our experiments, which include the Cox-Ingersoll-
Ross model and additive noise SDEs (noisy anharmonic oscillator, stochastic
FitzHugh-Nagumo model, underdamped Langevin dynamics), the new splitting
methods exhibit convergence rates of $O(h^{3/2})$ and outperform schemes
previously proposed in the literature.

###### keywords:

Numerical methods for SDEs, high order strong convergence, operator splitting

60H35, 60J65, 60L90, 65C30

##  1 Introduction

Stochastic differential equations (SDEs) are commonly used for modelling
random continuous -time phenomena, with applications ranging from finance [7,
57] and statistical physics [50, 59] to machine learning [44, 45, 75, 79, 82].
In such applications, SDE solutions can rarely be obtained exactly or in
closed-form, and so numerical methods and Monte Carlo simulation are often
employed in practice.

Figure 1: In the Monte Carlo paradigm, information about the Brownian motion
is generated and then mapped to a numerical solution of the SDE. Typically,
only Brownian increments are sampled.

In this paper, we present a study of high order splitting-based numerical
methods for Stratonovich SDEs of the form

(1) $dy_{t}=f(y_{t})\hskip 0.7113ptdt+g(y_{t})\circ dW_{t}\hskip
0.7113pt,\hskip 7.11317pty_{0}\in L^{2}(\mathbb{R}^{e}),$

where $W=(W^{1},\cdots,W^{d})=\\{W_{t}\\}_{t\in[0,T]}$ denotes a
$d$-dimensional Brownian motion, $L^{2}(\mathbb{R}^{e})$ is the space of
$\mathbb{R}^{e}$-valued square-integrable random variables, the vector fields
are given by $f\in\mathcal{C}^{2}(\mathbb{R}^{e},\mathbb{R}^{e})$ and
$g=(g_{1},\cdots,g_{d})\in\mathcal{C}^{3}(\mathbb{R}^{e},\mathbb{R}^{e\times
d})$ where we understand $g(y_{t})\circ dW_{t}=\sum_{i=1}^{d}g_{i}(y_{t})\circ
dW_{t}^{i}$. The columns $\\{g_{i}\\}_{1\leq i\leq d}$ of $g$ can each be
viewed as a vector field on $\mathbb{R}^{e}$ and are assumed to satisfy the
following commutativity condition:

(2) $\hskip 7.11317ptg_{i}^{\hskip 0.7113pt\prime}(y)g_{j}(y)=g_{j}^{\hskip
0.7113pt\prime}(y)g_{i}(y),\hskip 7.11317pt\forall y\in\mathbb{R}^{e}.$

We also assume $g_{i}$ are globally Lipschitz continuous with globally
Lipschitz derivatives.

Without the condition (2), high order numerical methods for SDEs require the
use, or approximation, of second iterated integrals of the Brownian motion
[71]. Generating both the increments and iterated integrals, or equivalently
Lévy areas, of Brownian motion is a difficult problem [16, 23] and beyond the
scope of this paper. We refer the reader to [20, 26, 27, 33, 64, 80] for
studies on Lévy area approximation. Nevertheless, there is a large variety of
SDEs used in applications that satisfy (2), such as SDEs with scalar, diagonal
or additive noise types. While we focus on schemes for SDEs satisfying the
commutativity condition (2), the error analysis that we introduce for
establishing convergence is generic and does not rely on this condition.

Inspired by rough path theory [32], which views SDEs as functions that map
Brownian motion to continuous paths (see Figure 1), we propose an
approximation $y^{\gamma}=\\{\hskip 0.7113pty_{r}^{\gamma}\hskip
0.7113pt\\}_{r\in[0,1]}\hskip 0.7113pt$ for (1) that comes from the controlled
differential equation (CDE),

(3) $\hskip
56.9055ptdy^{\gamma}_{r}=f(y^{\gamma}_{r})\,d\gamma^{\tau}(r)+g(y^{\gamma}_{r})\,d\gamma^{\hskip
0.7113pt\omega}(r)\hskip 0.7113pt,\hskip 28.45274pty^{\gamma}_{0}=y_{0}\hskip
0.7113pt,$

or equivalently

$y^{\gamma}_{r}=y^{\gamma}_{0}+\int_{0}^{r}f(y^{\gamma}_{u})\,d\gamma^{\tau}(u)+\int_{0}^{r}g(y^{\gamma}_{u})\,d\gamma^{\hskip
0.7113pt\omega}(u)\hskip 0.7113pt,$

where $\gamma=(\gamma^{\tau},\gamma^{\hskip
0.7113pt\omega})^{\top}:[0,1]\rightarrow\mathbb{R}^{1+d}$ is a parameterised
(continuous) piecewise linear path designed to match certain iterated
integrals of the “space-time” Brownian motion
$\\{(t,W_{t})\\}_{t\in[0,T]}\hskip 0.7113pt$. Since the path $\gamma$ is
piecewise linear, it immediately follows that
$d\gamma(r)=\frac{1}{r_{i+1}-r_{i}}\hskip 0.7113pt\gamma_{r_{i},r_{i+1}}\hskip
0.7113ptdr$ for $r\in[r_{i},r_{i+1}]$, where $r_{i}\in[0,1]$ is the parameter
value at the start of the $i$-th piece of $\gamma$ and
$\gamma_{r_{i},r_{i+1}}$ is the increment of the linear piece. Therefore the
CDE (3) reduces to a sequence of ODEs, corresponding to each piece of
$\gamma$, which can be discretised by a suitable ODE solver, such a Runge-
Kutta method. Furthermore, we will see that this approach can be interpreted
as a splitting method. We refer the reader to Section 3 of [8] for an overview
of splitting methods for SDEs.

More generally, CDEs are one of the key objects in rough path theory [31, 32,
55] (often referred to as “rough” differential equations). However, we
emphasise that this manuscript is not a rough paths paper – and no $p\hskip
0.7113pt$-variation or lift maps are used. Instead, we will heavily draw upon
ideas and interpretations from rough path theory. Similarly, we point towards
[11, 29, 43, 46, 56, 62, 63, 68] as works presenting results for stochastic
processes or continuous data streams, without “rough path” statements, but
making use of the machinery and insights that are provided by rough path
theory.

Perhaps the simplest example of an approximation with the form of the CDE (3)
is the Wong-Zakai approximation [68, 74, 81] where $\gamma$ is the standard
piecewise linear discretisation of space-time Brownian motion. However, the
Wong-Zakai approach only uses the increments of the Brownian motion, and is
thus constrained to a first order convergence rate for SDEs satisfying the
commutativity condition (2) [16]. We will show that by generating both
increments and “space-time” Lévy areas of the Brownian path, we can construct
paths $\gamma$ yielding order 3/2 strong convergence rates.

######  Definition 1.1.

The rescaled space-time Lévy area of a Brownian motion $W$ over an interval
$[s,t]$ corresponds to the signed area of the associated bridge process.

$\displaystyle
H_{s,t}:=\frac{1}{h}\int_{s}^{t}\Big{(}W_{s,u}-\frac{u-s}{h}\,W_{s,t}\Big{)}\,du\hskip
0.7113pt,$

where $h:=t-s$ and $W_{s,u}:=W_{u}-W_{s}$ for $u\in[s,t]$.

Figure 2: Space-time Lévy area gives the area between a Brownian path and its
linear approximant.

######  Remark 1.2.

It was shown in [29] that
$H_{s,t}\sim\mathcal{N}\big{(}0,\frac{1}{12}h\big{)}$ is independent of
$W_{s,t}$ when $d=1$. Since the coordinate processes of a Brownian motion are
independent, it therefore follows that $\hskip 0.7113ptW_{s,t}\hskip
0.7113pt\sim\hskip 0.7113pt\mathcal{N}\big{(}0,hI_{d}\big{)}\hskip 0.7113pt$
and $\hskip 0.7113ptH_{s,t}\hskip 0.7113pt\sim\hskip
0.7113pt\mathcal{N}\big{(}0,\frac{1}{12}hI_{d}\big{)}\hskip 0.7113pt$ are
independent.

Using the increment $W_{s,t}$ and space-time Lévy area $H_{s,t}$ of the
Brownian motion, we give an example of a path $\gamma$ and its associated 3/2
strong order spitting method.

######  Example 1.3.

Let $\gamma=(\gamma^{\tau},\gamma^{\hskip
0.7113pt\omega})^{\top}:[0,1]\rightarrow\mathbb{R}^{1+d}$ denote a piecewise
linear path where the vertices between pieces are at $\hskip
0.7113ptr_{i}:=\frac{i}{5}\hskip 0.7113pt$ for $\hskip 0.7113pt0\hskip
0.7113pt\leq\hskip 0.7113pti\hskip 0.7113pt\leq\hskip 0.7113pt5\hskip
0.7113pt$ and the increments are

$\displaystyle\gamma_{r_{i},r_{i+1}}=\begin{cases}\big{(}\frac{3-\sqrt{3}}{6}h,0\big{)},&\text{if
}\,i=0\\\\[6.0pt] \big{(}0,\frac{1}{2}W_{s,t}+\sqrt{3}H_{s,t}\big{)},&\text{if
}\,i=1\\\\[6.0pt] \big{(}\frac{\sqrt{3}}{3}h,0\big{)},&\text{if
}\,i=2\\\\[6.0pt]
\big{(}0,\frac{1}{2}W_{s,t}-\sqrt{3}H_{s,t}\big{)},&\textrm{if
}\,i=3\\\\[6.0pt] \big{(}\frac{3-\sqrt{3}}{6}h,0\big{)},&\text{if
}\,i=4.\end{cases}\begin{matrix}\\\ \hskip
14.22636pt\includegraphics[width=195.12767pt]{"Splitting_example".pdf}\end{matrix}$

Therefore, by replacing the driving signal $t\mapsto(t,W_{t})$ in the SDE (1)
with the parameterisation $r\mapsto\gamma_{r}$, the approximating CDE (3)
reduces to the splitting (formulated in a more classical way):

$\displaystyle y_{1}^{\gamma}$
$\displaystyle=\exp\bigg{(}\frac{3-\sqrt{3}}{6}\hskip 0.7113ptf(\hskip
0.7113pt\cdot\hskip 0.7113pt)\hskip 0.7113pth\bigg{)}\exp\bigg{(}g(\hskip
0.7113pt\cdot\hskip
0.7113pt)\Big{(}\frac{1}{2}W_{s,t}-\sqrt{3}H_{s,t}\Big{)}\bigg{)}$
$\displaystyle\hskip 28.45274pt\exp\bigg{(}\frac{\sqrt{3}}{3}\hskip
0.7113ptf(\hskip 0.7113pt\cdot\hskip 0.7113pt)\hskip
0.7113pth\bigg{)}\exp\bigg{(}g(\hskip 0.7113pt\cdot\hskip
0.7113pt)\Big{(}\frac{1}{2}W_{s,t}+\sqrt{3}H_{s,t}\Big{)}\bigg{)}\exp\bigg{(}\frac{3-\sqrt{3}}{6}\hskip
0.7113ptf(\hskip 0.7113pt\cdot\hskip 0.7113pt)\hskip
0.7113pth\bigg{)}y_{0}^{\gamma}\hskip 0.7113pt,$

where $\exp(V)\hskip 0.7113ptx$ denotes the solution $z_{1}$ at time $u=1$ of
the ODE $\,z^{\prime}=V(z),\,\,z(0)=x$.

The main result of this paper will allow us to establish the high order strong
convergence of (3) primarily by checking that the path $\gamma$ has the
following properties:

(4) $\displaystyle\hskip 5.69054pt\gamma^{\hskip
0.7113pt\omega}(1)-\gamma^{\hskip 0.7113pt\omega}(0)=W_{s,t}\hskip
0.7113pt,\hskip 21.33955pt\int_{0}^{1}\big{(}\gamma^{\hskip
0.7113pt\omega}(r)-\gamma^{\hskip 0.7113pt\omega}(0)\big{)}\hskip
0.7113ptd\gamma^{\tau}(r)=\int_{s}^{t}W_{s,u}\,du,$ (5)
$\displaystyle\mathbb{E}\bigg{[}\int_{0}^{1}\big{(}\gamma^{\hskip
0.7113pt\omega}(r)-\gamma^{\hskip 0.7113pt\omega}(0)\big{)}^{\otimes
2}d\gamma^{\tau}(r)\bigg{]}=\mathbb{E}\bigg{[}\int_{s}^{t}W_{s,u}^{\otimes
2}\,du\bigg{]}=\frac{1}{2}\hskip 0.7113pth^{2}I_{d}\hskip 0.7113pt.\hskip
21.33955pt$

These quantities correspond to terms in the Taylor expansions of the SDE (1)
and its CDE approximation (3). Provided $\gamma$ satisfies mild regularity
conditions, ensuring that the integrals in (4) and expected integrals in (5)
coincide allows the Taylor expansions (and hence solutions) of (1) and (3) to
be close in an $L^{2}(\mathbb{P})$ sense. Moreover, we will see that the
conditions (4) and (5) imposed on $\gamma$ are a specific instance of a more
general framework, which we shall make rigorous in Theorem 3.10. In terms of
proof methodologies for the error analysis, it is worth noting that ours
differs from previous works on splitting methods for SDEs [8, 28, 41, 50, 60,
66, 67, 78], which are either extensions of the Strang splitting [77] or use
the Baker-Campbell-Hausdorff formula for expanding the compositions of ODEs
(see [60] for the latter). That said, our approach does draw inspiration from
“Cubature on Wiener Space” [56] where (deterministic) piecewise linear paths
are used to weakly approximate SDEs. For some perspective, we present an
informal version of our main result, Theorem 3.10, which describes our
approach to high order splitting methods for commutative SDEs.

######  Theorem 1.4 (Convergence of path-based splitting for SDEs (informal
version)).

Given a fixed number of steps $N$, we will define a numerical solution
$Y=\\{Y_{k}\\}_{0\hskip 0.7113pt\leq\hskip 0.7113ptk\hskip 0.7113pt\leq N}$
for the SDE (1) over the finite time horizon $[0,T]$ as follows,

$Y_{k+1}:=\big{(}\text{Solution at time }r=1\text{ of CDE (\ref{eq:intro_CDE})
driven by }\gamma_{k}:[0,1]\rightarrow\mathbb{R}^{1+d}\hskip
0.7113pt\big{)}\big{(}Y_{k}\big{)}\hskip 0.7113pt,$

where each piecewise linear path $\gamma_{k}$ is constructed from
$\big{\\{}W_{t}:t\in\big{[}\frac{kT}{N},\frac{(k+1)T}{N}\big{]}\big{\\}}$, is
sufficiently regular (see Assumption A1), and for some fixed
$p\in\\{\frac{m}{2}\\}_{m\in\mathbb{N}}\hskip 0.7113pt$ satisfies
1\. the iterated integrals of $\gamma_{k}$ and $(t,W_{t})$ with order less
than $p-\frac{1}{2}$ coincide,
2\. the iterated integrals of $\gamma_{k}$ and $(t,W_{t})$ with order $p$
match in expectation.
Then, there exists a constant $C>0$, such that for sufficiently small
$h=\frac{T}{N}$, we have

(6) $\mathbb{E}\Big{[}\|Y_{k}-y_{kh}\|^{2}\Big{]}^{\frac{1}{2}}\leq
Ch^{p-\frac{1}{2}}\hskip 0.7113pt.$

for $k\in\\{1,\cdots,N\\}$. If $p=2$, and the SDE satisfies the commutativity
condition (2), then the estimate (6) holds under the assumption that each
$\gamma_{k}$ is sufficiently regular and has coordinate processes
$\\{\gamma_{k}^{\hskip 0.7113pt\omega,i}\\}_{1\hskip 0.7113pt\leq\hskip
0.7113pti\hskip 0.7113pt\leq d}$ that are independent, symmetric and satisfy

(7) $\displaystyle\gamma^{\hskip 0.7113pt\omega,i}_{k}(1)-\gamma^{\hskip
0.7113pt\omega,i}_{k}(0)=W_{kh,(k+1)h}^{i}\hskip 0.7113pt,\hskip
14.22636pt\gamma^{\tau}_{k}(1)-\gamma^{\tau}_{k}(0)=h,$ (8)
$\displaystyle\int_{0}^{1}\big{(}\gamma_{k}^{\hskip
0.7113pt\omega,i}(r)-\gamma_{k}^{\hskip 0.7113pt\omega,i}(0)\big{)}\hskip
0.7113ptd\gamma_{k}^{\tau}(r)=\int_{kh}^{(k+1)h}W_{kh,u}^{i}\,du\hskip
0.7113pt,$ (9)
$\displaystyle\mathbb{E}\bigg{[}\int_{0}^{1}\big{(}\gamma_{k}^{\hskip
0.7113pt\omega,i}(r)-\gamma_{k}^{\hskip
0.7113pt\omega,i}(0)\big{)}^{2}d\gamma_{k}^{\tau}(r)\bigg{]}=\frac{1}{2}\hskip
0.7113pth^{2}.\hskip 14.22636pt$

The paper is outlined as follows. In Section 2, we impose regularity
conditions on the path $\gamma$ and establish a (fourth) moment bound for the
solution of the CDE (3). In Section 3, we will present Taylor expansions for
both the SDE (1) and CDE (3). Since the error analysis of stochastic Taylor
approximations is well known, our focus is to obtain $L^{2}(\mathbb{P})$
estimates for the remainder terms in the CDE Taylor expansion. We end Section
3.3 with our main result, Theorem 3.10, which establishes convergence rates
for SDE splitting methods corresponding to the CDE (3) driven by a path
$\gamma$. To establish $3/2$ strong convergence rates in the setting where the
SDE satisfies the commutativity condition (2), we simplify certain terms in
the Taylor expansions of (1) and (3). This is the focus of Section 3.4, with
technical details in Appendix A.

In Section 4, we provide several examples of piecewise linear paths that
correspond to high order splitting methods. Two of the paths will be
constructed using a recently developed approximation [26, Theorem 5.1.2] for
the stochastic integral $\int_{s}^{t}W_{s,u}^{2}\hskip 0.7113ptdu$. To keep
this article self contained, we shall present the derivation of this new
integral approximation in Appendix B. Just like the approach of [29], this new
estimator can be obtained as a certain conditional expectation and is thus
optimal in an $L^{2}(\mathbb{P})$ sense. However, unlike in [29], we shall
generate an additional Rademacher random variable to improve the approximation
to aid in the construction of the piecewise linear paths.

In Section 5, we test the proposed splitting methods on some well-known SDEs,
including the Cox-Ingersoll-Ross [2, 17] and stochastic FitzHugh-Nagumo [8]
models. Due to the analytic tractability of these SDEs, our high order
splitting schemes will produce ODEs that can be either solved exactly or
further split into “solvable” ODEs. We also discuss the setting of additive
noise SDEs, where the drift vector field may not give analytically tractable
ODEs. For these problems, we propose applying a certain second order Runge-
Kutta method (Ralston’s method) to the “non-diffusion” ODE. Furthermore, we
show that the Taylor expansion of the resulting stochastic Ralston method
contains the high order terms needed to obtain order 3/2 strong convergence.
In the specific application to underdamped Langevin dynamics [50], we briefly
discuss the choice of splitting path and (third order) Runge-Kutta method
which can lead to a third order strong convergence rate to the SDE solution
(see [30] for further details). Finally, we demonstrate the improved accuracy
of the stochastic Ralston method when compared to the SRA1 scheme in [71], for
simulating a simple anharmonic oscillator. Moreover, as both methods require
two drift evaluations per step, we would expect stochastic Ralston to be
state-of-the-art for additive noise SDEs with expensive drifts. At the same
time, we also show that a simple splitting-based adjustment can improve the
accuracy of the standard Euler-Maruyama method for SDEs with additive noise.
We shall conclude the paper in Section 6 and discuss potential areas of future
research.

###  1.1 Notation

In this section, we summarise some of the notation in the paper. Given vectors
$a\in\mathbb{R}^{n}$ and $b\in\mathbb{R}^{m}$, we shall denote their tensor
product $a\otimes b\in\mathbb{R}^{nm}$ by $a\otimes b:=\\{a_{i}\hskip
0.7113ptb_{j}\\}_{i\hskip 0.7113pt=\hskip 0.7113pt1,\hskip 0.7113ptj\hskip
0.7113pt=\hskip 0.7113pt1}^{n,\hskip 10.38522ptm}$. We define iterated
integrals of $W$ and $\gamma:[0,1]\rightarrow\mathbb{R}^{n}$ as

(10) $\displaystyle I_{\alpha}(F)$
$\displaystyle:=\int_{0}^{h}\int_{0}^{r_{1}}\cdots\int_{0}^{r_{n-1}}F(y_{r_{n}})\,dB^{\alpha_{1}}_{r_{n}}\hskip
0.7113ptdB^{\alpha_{2}}_{r_{n-1}}\cdots
dB^{\alpha_{n-1}}_{r_{2}}dB^{\alpha_{n}}_{r_{1}},$ (11) $\displaystyle
I^{\gamma}_{\alpha}(F)$
$\displaystyle:=\int_{0}^{1}\int_{0}^{r_{1}}\cdots\int_{0}^{r_{n-1}}F(y^{\gamma}_{r_{n}})\,d\gamma^{\alpha_{1}}(r_{n})\hskip
0.7113pt\cdots d\gamma^{\alpha_{n}}(r_{1})\hskip 0.7113pt,$

where $y$ and $y^{\gamma}$ are the solutions of the SDE (1) and its CDE
approximation (3), $F:\mathbb{R}^{e}\rightarrow\mathbb{R}^{e}$,
$\alpha=(\alpha_{1},\cdots,\alpha_{n})\in\\{\tau,\omega\\}^{n}$ is a multi-
index, $dB^{\tau}_{r}=dr$, $dB^{\hskip 0.7113pt\omega}_{r}=\otimes\circ
dW_{r}\hskip 0.7113pt$ and, if $\gamma(r)=(\gamma^{\tau}(r),\gamma^{\hskip
0.7113pt\omega}(r))^{\top}$, then we write
$d\gamma(r)^{\tau}=d\gamma^{\tau}(r)$ and $d\gamma(r)^{\hskip
0.7113pt\omega}=\otimes\,d\gamma^{\hskip 0.7113pt\omega}(r)$. We denote the
set of multi-indices by $\mathcal{A}=\cup_{n\geq 0}\\{\tau,\omega\\}^{n}$. We
also define the integrals,

(12) $J_{\alpha}(F):=I_{\alpha}(F)-F(y_{0})I_{\alpha}(1)\hskip 0.7113pt,\hskip
14.22636ptJ^{\gamma}_{\alpha}(F):=I^{\gamma}_{\alpha}(F)-F(y^{\gamma}_{0})I^{\gamma}_{\alpha}(1)\hskip
0.7113pt,$

where we understand $I_{\alpha}(1)$ and $I^{\gamma}_{\alpha}(1)$ as defined in
(10) and (11), but with $F(y)$ replaced by the scalar $1$. For a given multi-
index $\alpha=(\alpha_{i})_{1\hskip 0.7113pt\leq\hskip 0.7113pti\hskip
0.7113pt\leq\hskip 0.7113ptn}\hskip 0.7113pt$, we will define its order by
$\textrm{ord}(\alpha):=|\alpha|_{\tau}+\frac{1}{2}|\alpha|_{\omega}\hskip
0.7113pt$, where $|\alpha|_{\tau}:=\sum_{i=1}^{n}\mathbf{1}_{\alpha_{i}\hskip
0.7113pt=\hskip 0.7113pt\tau}$ and
$|\alpha|_{\omega}:=\sum_{i=1}^{n}\mathbf{1}_{\alpha_{i}\hskip 0.7113pt=\hskip
0.7113pt\omega}\hskip 0.7113pt$.

Given normed vector spaces $U$ and $V$, $\mathcal{C}_{\mathrm{Lip}}^{p}(U,V)$
will be the subspace of $\mathcal{C}^{p}(U,V)$ containing globally Lipschitz
continuous functions with their $p$ derivatives globally Lipschitz continuous.
We write the Lipscthiz constant of a function $F$ as
$\|F\|_{\mathrm{Lip}\hskip 0.7113pt\text{-}1}$. The above notation shall be
employed in the study of Taylor expansions in Section 3.

Throughout, $\|\cdot\|$ will denote the standard Euclidean norm on
$\mathbb{R}^{n}$ and $L^{p}(\mathbb{R}^{n})$ is the space of
$\mathbb{R}^{n}$-valued random variables with finite $p\hskip 0.7113pt$-th
moments (i.e. $\mathbb{E}\big{[}\|X\|^{p}\big{]}<\infty$). Given a continuous
path $\gamma:[0,1]\rightarrow\mathbb{R}^{n}$, we will denote its length using
the notation,

$\|\gamma\|_{1\text{-var},[0,1]}\hskip 0.7113pt=\int_{0}^{1}|\hskip
0.7113ptd\gamma(r)|\hskip 0.7113pt:=\sup_{\begin{subarray}{c}0\hskip
0.7113pt=\hskip 0.7113ptr_{0}\hskip 0.7113pt<\hskip 0.7113ptr_{1}\hskip
0.7113pt<\hskip 0.7113pt\cdots\hskip 0.7113pt<\hskip 0.7113ptr_{N}\hskip
0.7113pt=\hskip 0.7113pt1,\\\\[3.0pt] N\hskip 0.7113pt\geq\hskip
0.7113pt0.\end{subarray}}\bigg{(}\sum_{i=0}^{N-1}\|\gamma(r_{i+1})-\gamma(r_{i})\|\bigg{)}\hskip
0.7113pt.$

When defining numerical methods, we sometimes use $W_{k}$ as shorthand for
$W_{t_{k}\hskip 0.7113pt,t_{k+1}}$, (and similarly $H_{k}$ and $n_{k}$ instead
of $H_{t_{k}\hskip 0.7113pt,t_{k+1}}$ and $n_{t_{k}\hskip 0.7113pt,t_{k+1}}$).

##  2 Main assumption and preliminary results

Before we prove the strong convergence of the CDE (3), we will first establish
a moment bound for its solution. As discussed previously, this requires us to
make certain assumptions on the path $\gamma$. Our main assumption (given
below) ensures the path $\gamma$ scales like Brownian motion.

######  Assumption A1 (Brownian-like scaling).

Let $\gamma=(\gamma^{\tau},\gamma^{\hskip
0.7113pt\omega})^{\top}:[0,1]\to\mathbb{R}^{1+d}$ be a piecewise linear path
with $m\in\mathbb{N}$ components that have, almost surely, finite length. For
$i\geq 0$, we denote the increment of the $i$-th piece of $\gamma$ by
$\gamma_{r_{i},r_{i+1}}$ and assume that

  1. 1.

$\gamma^{\tau}_{r_{i},r_{i+1}}\hskip 0.7113pt$, the increment in the time
component of $\gamma$, is deterministic.

  2. 2.

$\gamma^{\tau}_{r_{i},r_{i+1}}$ scales with the step size $h$ and the
increment in the space component, $\gamma^{\hskip
0.7113pt\omega}_{r_{i},r_{i+1}}$, has finite even moments scaling with $h$.
Concretely, we have

$\displaystyle\gamma^{\tau}_{r_{i},r_{i+1}}=O(h),\hskip
7.11317pt\textrm{and}\quad\mathbb{E}\big{[}|(\gamma^{\hskip
0.7113pt\omega}_{r_{i},r_{i+1}})_{j}|^{2k}\big{]}=O(h^{k}),$

for every $j\in\\{1,\cdots,d\\}$.

######  Remark 2.1 (Comment on Assumption A1).

We impose that $\gamma^{\tau}$ is deterministic for convenience and,
inspecting the proof, one may be able to lift this constraint. Moreover, we
expect our methodology can accommodate for randomised algorithms (see [6, 39,
48, 49, 73] for examples of SDE solvers with a randomised time component).

We now present the main result of this section – a moment bound for the CDE,
which will be used to control remainder terms of the Taylor expansion
discussed later. Following the approach of [32, Theorem 3.7], we obtain our
main result, Theorem 2.2.

######  Theorem 2.2 (Fourth moment bound for CDEs).

Let $\gamma$ satisfy Assumption A1 and let $y^{\gamma}$ denote the solution to
(3) with $y^{\gamma}_{0}\in L^{4}(\mathbb{R}^{e})$. Suppose that $f$ and $g$
satisfy

(13) $\displaystyle\|f(y)\|\leq C(1+\|y\|),\quad\textrm{and}\quad\|g(y)\|\leq
C(1+\|y\|),$

with $\mathbb{E}\big{[}\exp\big{(}16\hskip
0.7113ptC\int_{0}^{1}|d\gamma(u)|\big{)}\big{]}<\infty$. Then there exists a
positive constant $\widetilde{C}>0$, depending only the path $\gamma$ and the
growth constant $C$ in (13), such that for $r\in[0,1]$,

(14) $\displaystyle\mathbb{E}\big{[}\hskip
0.7113pt\|y_{r}^{\gamma}-y_{0}^{\gamma}\|^{4}\big{]}\leq\widetilde{C}\hskip
0.7113pth^{2}\big{(}1+\mathbb{E}\big{[}\|y_{0}\|^{4}\big{]}\big{)}\hskip
0.7113pt.$

######  Proof 2.3.

Let $G:\mathbb{R}^{e}\to\mathbb{R}^{e\times(d+1)}$ have first column given by
$f:\mathbb{R}^{e}\to\mathbb{R}^{e}$ and the rest of the matrix given by
$g:\mathbb{R}^{e}\to\mathbb{R}^{e\times d}$. Then the growth assumption (13)
implies that $\|G(y)\|\leq C(1+\|y\|)$. Thus, by direct application of [32,
Theorem 3.7], we have

$\displaystyle\|y^{\gamma}_{r}-y^{\gamma}_{0}\|\leq
C(1+\|y^{\gamma}_{0}\|)\exp\left(2\hskip
0.7113ptC\int_{0}^{r}|d\gamma(u)|\right)\int_{0}^{r}|d\gamma(u)|\hskip
0.7113pt,$

for $r\in[0,1]$. Consequently,

$\displaystyle\mathbb{E}\big{[}\|y^{\gamma}_{r}-y^{\gamma}_{0}\|^{4}\big{]}$
$\displaystyle\leq C^{4}\hskip
0.7113pt\mathbb{E}\left[\big{(}1+\|y^{\gamma}_{0}\|\big{)}^{4}\exp\left(8\hskip
0.7113ptC\int_{0}^{r}|d\gamma(u)|\right)\left(\int_{0}^{r}|d\gamma(u)|\right)^{4}\right]$
$\displaystyle\leq C^{4}\hskip
0.7113pt\big{(}1+\mathbb{E}\big{[}\|y_{0}^{\gamma}\|^{4}\big{]}\big{)}\mathbb{E}\bigg{[}\exp\bigg{(}16\hskip
0.7113ptC\int_{0}^{r}|d\gamma(u)|\bigg{)}\bigg{]}^{\frac{1}{2}}\mathbb{E}\Bigg{[}\bigg{(}\int_{0}^{r}|d\gamma(u)|\bigg{)}^{8}\Bigg{]}^{\frac{1}{2}}$
$\displaystyle\leq\widetilde{C}\hskip
0.7113pth^{2}\big{(}1+\mathbb{E}\big{[}\|y_{0}\|^{4}\big{]}\big{)},$

where we used the independence of $y_{0}$ and $\gamma$, the Cauchy-Schwarz
inequality and

$\displaystyle\mathbb{E}\Bigg{[}\bigg{(}\int_{0}^{r}|d\gamma(u)|\bigg{)}^{8}\hskip
0.7113pt\Bigg{]}\leq\mathbb{E}\Bigg{[}\bigg{(}\int_{0}^{1}|d\gamma(u)|\bigg{)}^{8}\hskip
0.7113pt\Bigg{]}$
$\displaystyle=\mathbb{E}\Bigg{[}\bigg{(}\sum_{i=0}^{m-1}\big{\|}\gamma_{r_{i}\hskip
0.7113pt,r_{i+1}}\big{\|}\bigg{)}^{8}\hskip 0.7113pt\Bigg{]}$
$\displaystyle\leq m^{7}\sum_{i=0}^{m-1}\mathbb{E}\Big{[}\|\gamma_{r_{i}\hskip
0.7113pt,r_{i+1}}\big{\|}^{8}\Big{]}=O(h^{4}),$

which follows by Jensen’s inequality and Assumption A1.

##  3 Taylor expansions and error analysis

We consider the Taylor expansions of both the Stratonovich SDE (1) and the CDE
(3) driven by a splitting path $\gamma$. By matching the lower order terms in
the Taylor expansions and showing that the remainder terms are higher order,
we can bound local errors for our splitting schemes. We then apply Milstein
and Tretyakov’s framework for mean-square error analysis [59] to obtain a
global strong convergence rate – which is our main result in Theorem 3.10.

###  3.1 Stratonovich Taylor expansion

Letting $y$ denote the solution to (1), we have the following chain rule (see
[32, Exercise 3.17]) for any function $F\in\mathcal{C}^{1}(\mathbb{R}^{e})$,

(15) $F(y_{r})=F(y_{0})+\int_{0}^{r}F^{\prime}(y_{s})\circ dy_{s}\hskip
0.7113pt,$

By expanding “$dy_{s}$” and iteratively applying (15), we obtain the Taylor
expansion.

######  Proposition 3.1 (Stochastic Taylor expansion of the Stratonovich SDE
(1) [5, Proposition 1.1], [47, Theorem 5.6.1]).

Let $p\in\\{\frac{k}{2}\\}_{k\in\mathbb{N}}$,
$f\in\mathcal{C}_{\mathrm{Lip}}^{\lceil
p-1\rceil}(\mathbb{R}^{e},\mathbb{R}^{e})$ and
$g\in\mathcal{C}_{\mathrm{Lip}}^{2p-1}(\mathbb{R}^{e},\mathbb{R}^{e\times
d})$. The Stratonovich Taylor expansion of (1), up to order $p$, is

(16) $\displaystyle y_{h}=y_{0}+\sum_{\begin{subarray}{c}\alpha\hskip
0.7113pt\in\mathcal{A}\,,\\\\[2.0pt]
\textrm{ord}(\alpha)\leq\,p\end{subarray}}V(\alpha)(y_{0})I_{\alpha}(1)+R_{p}(h,y_{0}),$

where, we recall the definition of
$\textrm{ord}(\alpha):=|\alpha|_{\tau}+\frac{1}{2}|\alpha|_{\omega}\hskip
0.7113pt$ after equation (12), and

(17) $\displaystyle R_{p}(h,y_{0}):=\sum_{\begin{subarray}{c}\alpha\hskip
0.7113pt\in\mathcal{A}\,,\\\\[2.0pt]
\textrm{ord}(\alpha)=p\end{subarray}}J_{\alpha}(V(\alpha))\hskip 0.7113pt,$

with the vector field derivatives $V(\alpha):\mathbb{R}^{e}\to
L((\mathbb{R}^{d})^{\otimes|\alpha|_{\omega}},\mathbb{R}^{e})$ defined for
multi-indices recursively by $V(\tau)(y):=f(y)$, $V(\omega)(y):=g(y)$ and

$V(l\beta)(y)=V(\beta)^{\prime}V(l)(y),$

where $l\in\\{\tau,\omega\\}$ and $l\beta:=(l,\beta_{1},\cdots,\beta_{n})$
denotes concatenation. Moreover, we have

(18)
$\mathbb{E}\big{[}\|R_{p}(h,y_{0})\|^{2}\big{]}^{1/2}=O(h^{p+\frac{1}{2}})\hskip
0.7113pt.$

###  3.2 Controlled Taylor expansion

We now present a CDE Taylor expansion. Just as with the Stratonovich SDE, we
have the following chain rule for $F\in\mathcal{C}^{1}(R^{e})$,

(19) $\hskip
14.22636ptF(y^{\gamma}_{r})=F(y^{\gamma}_{0})+\int_{0}^{r}F^{\prime}(y^{\gamma}_{s})\,dy^{\gamma}_{s}\hskip
0.7113pt,\hskip 14.22636ptr\in[0,1]\hskip 0.7113pt,$

where $y^{\gamma}$ denotes the solution to the CDE (3). Again, just as in the
SDE setting, by expanding “$\hskip 0.7113ptdy^{\gamma}_{s}\hskip 0.7113pt$”
and iteratively applying (19), we can obtain a Taylor expansion.

######  Proposition 3.2.

Let $p\in\\{\frac{k}{2}\\}_{k\in\mathbb{N}}$,
$f\in\mathcal{C}_{\mathrm{Lip}}^{\lceil
p-1\rceil}(\mathbb{R}^{e},\mathbb{R}^{e})$ and
$g\in\mathcal{C}_{\mathrm{Lip}}^{2p-1}(\mathbb{R}^{e},\mathbb{R}^{e\times
d})$. Then the (controlled) Taylor expansion of the CDE (3) up to order $p$ is
given by

(20) $y^{\gamma}_{1}=y^{\gamma}_{0}+\sum_{\begin{subarray}{c}\alpha\hskip
0.7113pt\in\mathcal{A}\hskip 0.7113pt,\\\\[2.0pt] \textrm{ord}(\alpha)\hskip
0.7113pt\leq\,p\end{subarray}}V(\alpha)(y^{\gamma}_{0})I^{\gamma}_{\alpha}(1)+R^{\gamma}_{p}(h,y_{0}^{\gamma})\hskip
0.7113pt,$

where, using the same notation for vector field derivatives as Proposition
3.1, we have

(21) $R^{\gamma}_{p}(h,y_{0}^{\gamma}):=\sum_{\begin{subarray}{c}\alpha\hskip
0.7113pt\in\mathcal{A}\hskip 0.7113pt,\\\\[2.0pt] \textrm{ord}(\alpha)\hskip
0.7113pt=\,p\end{subarray}}J^{\gamma}_{\alpha}(V(\alpha))\hskip 0.7113pt.$

######  Remark 3.3.

We note the proof of Proposition 3.2 is essentially identical to that of the
Stratonovich Taylor expansion (but with $t\mapsto(t,W_{t})$ replaced by
$r\mapsto(\gamma_{r}^{\tau},\gamma_{r}^{\hskip 0.7113pt\omega})$). Moreover,
we consider SDEs defined by Stratonovich integration precisely because it
results in the “same” integration by parts formula as for Riemann-Stieltjes
integrals.

We now show that the size of the terms in the CDE expansion (20) are dictated
by the order of $\alpha$. This will allow us to obtain a bound on the
remainder term $R^{\gamma}_{p}$ which we will then use to establish the
(strong) convergence rate of the CDE approximation. To begin, we shall
consider the iterated integrals $I^{\gamma}_{\alpha}\hskip 0.7113pt$, which do
not depend on $f$ or $g$.

######  Lemma 3.4.

Suppose the path $\gamma$ satisfies Assumption A1 and let
$\alpha\in\mathcal{A}$. Then

$\displaystyle\mathbb{E}\big{[}\|I^{\gamma}_{\alpha}(1)\|^{2}\big{]}=O\big{(}h^{2\hskip
0.7113pt\textrm{ord}(\alpha)}\big{)}\hskip 0.7113pt.$

######  Proof 3.5.

Since $\gamma$ is piecewise linear, we may split the iterated integral of
$\gamma$ into a finite sum of iterated integrals over intervals where the
derivatives of $\gamma$ are constant. This directly follows by the standard
additive property of Riemann-Stieltjes integrals. We may convert these path
integrals into regular (deterministic) integrals over these intervals as
$d\gamma^{\tau}(r)=\frac{1}{r_{i+1}-r_{i}}\gamma^{\tau}_{r_{i},r_{i+1}}dr$ and
$d\gamma^{\hskip 0.7113pt\omega}(r)=\frac{1}{r_{i+1}-r_{i}}\gamma^{\hskip
0.7113pt\omega}_{r_{i},r_{i+1}}dr$ for $r\in[r_{i},r_{i+1}]$. By Assumption
A1, $\gamma^{\tau}_{r_{i},r_{i+1}}=O(h)$ is a deterministic constant and we
therefore have

$\displaystyle\mathbb{E}\big{[}\|I^{\gamma}_{\alpha}(1)\|^{2}\big{]}$
$\displaystyle\leq
Ch^{2|\alpha|_{\tau}}\sum_{\textrm{intervals}}\mathbb{E}\bigg{[}\Big{\|}\bigotimes_{j=1}^{{|\alpha|_{\omega}}}\gamma^{\hskip
0.7113pt\omega}_{r^{j}_{i},r^{j}_{i+1}}\Big{\|}^{2}\bigg{]}\hskip 0.7113pt,$
$\displaystyle=Ch^{2|\alpha|_{\tau}}\sum_{\textrm{intervals}}\mathbb{E}\bigg{[}\sum_{i_{1}=1}^{d}\dots\sum^{d}_{i_{|\alpha|_{\omega}}=1}\big{(}\gamma^{\hskip
0.7113pt\omega}_{r^{1}_{i},r^{1}_{i+1}}\big{)}_{i_{1}}^{2}\times\cdots\times\big{(}\gamma^{\hskip
0.7113pt\omega}_{r^{|\alpha|_{\omega}}_{i},r^{|\alpha|_{\omega}}_{i+1}}\big{)}_{i_{|\alpha|_{\omega}}}^{2}\bigg{]}\hskip
0.7113pt,$

where “intervals” refers to the finite collection of subdomains of the simplex
with
$d\gamma^{\tau}(r)=\frac{1}{r_{i+1}-r_{i}}\gamma^{\tau}_{r_{i},r_{i+1}}dr$ and
$d\gamma^{\hskip 0.7113pt\omega}(r)=\frac{1}{r_{i+1}-r_{i}}\gamma^{\hskip
0.7113pt\omega}_{r_{i},r_{i+1}}dr$. We can then estimate the $\gamma^{\hskip
0.7113pt\omega}_{r_{i},r_{i+1}}$ terms by iteratively applying Hölder’s
inequality to the expectation and applying the assumption that
$\mathbb{E}\big{[}|(\gamma^{\hskip
0.7113pt\omega}_{r_{i},r_{i+1}})_{j}|^{2k}\big{]}=O(h^{k})$ for
$k\in\mathbb{N}$. This implies that

$\displaystyle\mathbb{E}\big{[}\|I^{\gamma}_{\alpha}(1)\|^{2}\big{]}\leq
C_{d,m,|\alpha|}h^{2\hskip 0.7113pt\textrm{ord}(\alpha)}.$

We now consider the $J^{\gamma}_{\alpha}$ terms, which will follow in much the
same way as for $I^{\gamma}_{\alpha}$.

######  Lemma 3.6.

Suppose that the assumptions of Theorem 2.2 hold. Let $\alpha\in\mathcal{A}$
and $F:\mathbb{R}^{e}\to L\big{(}(\mathbb{R}^{d})^{\otimes
b},\mathbb{R}^{e}\big{)}$ be a globally Lipschitz continuous map for some
$b\geq 1$, then

$\displaystyle\mathbb{E}\big{[}\|J^{\gamma}_{\alpha}(F)\|^{2}\big{]}=O\big{(}h^{2\hskip
0.7113pt\textrm{ord}(\alpha)+1}\big{)}\hskip 0.7113pt.$

######  Proof 3.7.

Just as in the previous proof, we may split the iterated integral of the
piecewise linear path $\gamma$ into a finite sum of iterated integrals over
the intervals where both
$d\gamma^{\tau}(r)=\frac{1}{r_{i+1}-r_{i}}\gamma^{\tau}_{r_{i},r_{i+1}}dr$ and
$d\gamma^{\hskip 0.7113pt\omega}(r)=\frac{1}{r_{i+1}-r_{i}}\gamma^{\hskip
0.7113pt\omega}_{r_{i},r_{i+1}}dr$ for all $r\in[r_{i},r_{i+1}]$. Applying
Jensen’s and Hölder’s inequalities to the finite sum in
$\mathbb{E}\big{[}\|J^{\gamma}_{\alpha}(F)\|^{2}\big{]}$ yields

$\displaystyle\mathbb{E}\big{[}\|J^{\gamma}_{\alpha}(F)\|^{2}\big{]}$
$\displaystyle\leq
Ch^{2|\alpha|_{\tau}}\sum_{\textrm{intervals}}\mathbb{E}\Bigg{[}\,\bigg{\|}\underset{0\hskip
0.7113pt<\hskip 0.7113ptr_{|\alpha|}\hskip 0.7113pt<\hskip
0.7113pt\cdots\hskip 0.7113pt<\hskip 0.7113ptr_{1}\hskip 0.7113pt<\hskip
0.7113pt1}{\idotsint}F(y^{\gamma}_{r_{|\alpha|}})-F(y^{\gamma}_{0})\hskip
0.7113ptdr_{|\alpha|}\cdots dr_{1}\bigg{\|}^{4}\Bigg{]}^{\frac{1}{2}}$
$\displaystyle\hskip
99.58464pt\times\mathbb{E}\Bigg{[}\bigg{\|}\bigotimes_{j=1}^{|\alpha|_{\omega}}\gamma^{\hskip
0.7113pt\omega}_{r^{j}_{i},r^{j}_{i+1}}\bigg{\|}^{4}\Bigg{]}^{\frac{1}{2}}.$

By applying Jensen’s inequality to the uniform distribution on $[s,t]$, we
have that $\big{\|}\int_{s}^{t}z_{r}\hskip
0.7113ptdr\big{\|}^{4}\leq(t-s)^{3}\int_{s}^{t}\|z_{r}\|^{4}\,dr$ for any
continuous integrable process $z_{r}\hskip 0.7113pt$. Therefore,

$\displaystyle\mathbb{E}\big{[}\|J^{\gamma}_{\alpha}(F)\|^{2}\big{]}$
$\displaystyle\leq
C_{1}h^{2|\alpha|_{\tau}}\sum_{\textrm{intervals}}\mathbb{E}\Bigg{[}\,\,\underset{0\hskip
0.7113pt<\hskip 0.7113ptr_{|\alpha|}\hskip 0.7113pt<\hskip
0.7113pt\cdots\hskip 0.7113pt<\hskip 0.7113ptr_{1}\hskip 0.7113pt<\hskip
0.7113pt1}{\idotsint}\big{\|}F(y^{\gamma}_{r_{|\alpha|}})-F(y^{\gamma}_{0})\big{\|}^{4}\hskip
0.7113ptdr_{|\alpha|}\cdots dr_{1}\Bigg{]}^{\frac{1}{2}}$ $\displaystyle\hskip
99.58464pt\times\mathbb{E}\Bigg{[}\bigg{\|}\bigotimes_{j=1}^{|\alpha|_{\omega}}\gamma^{\hskip
0.7113pt\omega}_{r^{j}_{i},r^{j}_{i+1}}\bigg{\|}^{4}\Bigg{]}^{\frac{1}{2}}.$

By repeatedly applying Hölder’s inequality, we can estimate the last term as
$O(h^{|\alpha|_{\omega}})$. Since
$\textrm{ord}(\alpha)=|\alpha|_{\tau}+\frac{1}{2}|\alpha|_{\omega}$, it
follows from the global Lipschitz continuity of $F$ that

$\displaystyle\mathbb{E}\big{[}\|J^{\gamma}_{\alpha}(F)\|^{2}\big{]}$
$\displaystyle\leq C_{2}\|F\|^{2}_{\mathrm{Lip}\text{-}1}h^{2\hskip
0.7113pt\textrm{ord}(\alpha)}$ $\displaystyle\hskip
28.45274pt\times\sum_{\textrm{intervals}}\Bigg{(}\,\,\underset{0\hskip
0.7113pt<\hskip 0.7113ptr_{|\alpha|}\hskip 0.7113pt<\hskip
0.7113pt\cdots\hskip 0.7113pt<\hskip 0.7113ptr_{1}\hskip 0.7113pt<\hskip
0.7113pt1}{\idotsint}\mathbb{E}\Big{[}\big{\|}y^{\gamma}_{r_{|\alpha|}}-y^{\gamma}_{0}\big{\|}^{4}\Big{]}\hskip
0.7113ptdr_{|\alpha|}\cdots dr_{1}\Bigg{)}^{\frac{1}{2}}.$

By Theorem 2.2, we have
$\mathbb{E}\big{[}\big{\|}y^{\gamma}_{r_{|\alpha|}}-y^{\gamma}_{0}\big{\|}^{4}\big{]}=O(h^{2})$
and thus the result follows.

###  3.3 Main result

Now that we have Taylor expansions for both the CDE and the Stratonovich SDE,
along with control over the size of the remainder terms in each, we can
establish the strong convergence properties of path-based splitting schemes.
We first obtain local strong and weak error estimates using a direct
application of Lemmas 3.4 and 3.6 before applying the framework of Milstein
and Tretyakov [59], which allows us to prove a global strong convergence rate
for the approximating CDE.

######  Theorem 3.8 (Local error estimates).

Suppose that the path $\gamma:[0,1]\to\mathbb{R}^{1+d}$ satisfies Assumption
A1 and for a fixed $p\in\\{\frac{k}{2}\\}_{k\in\mathbb{N}}\hskip 0.7113pt$,
let $f\in\mathcal{C}_{\mathrm{Lip}}^{\lceil
p-1\rceil}(\mathbb{R}^{e},\mathbb{R}^{e})$ and
$g\in\mathcal{C}_{\mathrm{Lip}}^{2p-1}(\mathbb{R}^{e},\mathbb{R}^{e\times
d})$. Suppose also that the assumptions of Theorem 2.2 hold and the integrals
$I^{\gamma}_{\alpha}(1)$ and $I_{\alpha}(1)$ agree almost surely for
$\alpha\in\mathcal{A}$ with $\textrm{ord}(\alpha)\leq p-\frac{1}{2}$ and in
expectation for all $\alpha\in\mathcal{A}$ with $\textrm{ord}(\alpha)=p$. Let
$Y_{1}$ denote an approximation (e.g. using an ODE solver) of the CDE solution
$\\{y^{\gamma}_{r}\\}_{r\in[0,1]}$ driven by $\gamma$, such that
$y_{0}^{\gamma}=y_{0}$ and

$\mathbb{E}\big{[}\|y^{\gamma}_{1}-Y_{1}\|\big{]}^{\frac{1}{2}}=O(h^{p}),\hskip
14.22636pt\text{and}\hskip 14.22636pt\big{\|}\mathbb{E}[\hskip
0.7113pty^{\gamma}_{1}]-\mathbb{E}[Y_{1}]\big{\|}=O\big{(}h^{p+\frac{1}{2}}\big{)},$

where $y=\\{y_{t}\\}_{t\in[0,h]}$ is the solution of the SDE (1) and $h>0$ is
the step size. Then

$\mathbb{E}\big{[}\|y_{h}-Y_{1}\|\big{]}^{\frac{1}{2}}=O(h^{p}),\hskip
14.22636pt\text{and}\hskip
14.22636pt\big{\|}\mathbb{E}[y_{h}]-\mathbb{E}[Y_{1}]\big{\|}=O\big{(}h^{p+\frac{1}{2}}\big{)}.$

######  Proof 3.9.

We start by proving the local strong error. By the triangle inequality,

$\mathbb{E}\big{[}\|y_{h}-Y_{1}\|^{2}\big{]}^{\frac{1}{2}}\leq\mathbb{E}\big{[}\|y_{h}-y^{\gamma}_{1}\|^{2}\big{]}^{\frac{1}{2}}+\mathbb{E}\big{[}\|y^{\gamma}_{1}-Y_{1}\|^{2}\big{]}^{\frac{1}{2}}=\mathbb{E}\big{[}\|y_{h}-y^{\gamma}_{1}\|^{2}\big{]}^{\frac{1}{2}}+O(h^{p}),$

as the second term is the difference between the CDE solution and its
approximation.

Recall the remainder terms $R_{p}(h,y_{0})$ and $R^{\gamma}_{p}(h,y_{0})$ in
Propositions 3.1 and 3.2. Then, by another two applications of the triangle
inequality, it directly follows that

$\displaystyle\mathbb{E}\big{[}\|y_{h}-Y_{1}\|^{2}\big{]}^{\frac{1}{2}}$
$\displaystyle\leq\mathbb{E}\big{[}\|(y_{h}-R_{p}(h,y_{0}))-(y^{\gamma}_{1}-R^{\gamma}_{p}(h,y_{0}))\|^{2}\big{]}^{\frac{1}{2}}$
$\displaystyle\hskip
14.22636pt+\mathbb{E}\big{[}\|R_{p}(h,y_{0})\|^{2}\big{]}^{\frac{1}{2}}+\mathbb{E}\big{[}\|R^{\gamma}_{p}(h,y_{0})\|^{2}\big{]}^{\frac{1}{2}}+O(h^{p}),$

where the first term is simply the difference in the Taylor expansions, up to
order $p$, of the SDE solution $y_{h}$ and the CDE solution $y_{1}^{\gamma}$.
Therefore, by the assumption that all integrals of the form
$I^{\gamma}_{\alpha}(1)$ are matched almost surely for
$\textrm{ord}(\alpha)\leq p-\frac{1}{2}\hskip 0.7113pt$, we have

$\mathbb{E}\big{[}\|y_{h}-Y_{1}\|^{2}\big{]}^{\frac{1}{2}}\leq\mathbb{E}\big{[}\|R_{p}(h,y_{0})\|^{2}\big{]}^{\frac{1}{2}}+\mathbb{E}\big{[}\|R^{\gamma}_{p}(h,y_{0})\|^{2}\big{]}^{\frac{1}{2}}+O(h^{p}).$

By Proposition 3.1, the SDE remainder term will satisfy
$\mathbb{E}\big{[}\|R_{p}(h,y_{0})\|^{2}\big{]}^{\frac{1}{2}}=O(h^{p+\frac{1}{2}})$.
On the other hand, $R^{\gamma}_{p}(h,y_{0})$ is given by (21) and therefore,
by Lemma 3.6, we have

$\mathbb{E}\big{[}\|R^{\gamma}_{p}(h,y_{0})\|^{2}\big{]}^{\frac{1}{2}}=O\big{(}h^{p+\frac{1}{2}}\big{)}.$

This gives the desired result for the local strong error, that
$\mathbb{E}\big{[}\|y_{h}-Y_{1}\|^{2}\big{]}^{\frac{1}{2}}=O(h^{p})$.

We now turn our attention to the local weak error. Using the triangle
inequality and the same Taylor expansions as in the proof of local strong
error, it follows that

$\displaystyle\big{\|}\mathbb{E}[y_{h}]-\mathbb{E}[Y_{1}]\big{\|}$
$\displaystyle\leq\big{\|}\mathbb{E}\big{[}y_{h}-R_{p}(h,y_{0})\big{]}-\mathbb{E}\big{[}y^{\gamma}_{1}-R^{\gamma}_{p}(h,y_{0})\big{]}\big{\|}$
$\displaystyle\hskip
14.22636pt+\big{\|}\mathbb{E}[y^{\gamma}_{1}]-\mathbb{E}[Y_{1}]\big{\|}+\big{\|}\mathbb{E}[R_{p}(h,y_{0})]\big{\|}+\big{\|}\mathbb{E}[R^{\gamma}_{p}(h,y_{0})]\big{\|}.$

From our assumption, the $I^{\gamma}_{\alpha}(1)$ terms in the SDE and CDE
Taylor expansions are matched in expectation for $\textrm{ord}(\alpha)\leq p$
and, therefore, the first term disappears. Moreover, we assume
$\|\mathbb{E}[y_{h}]-\mathbb{E}[Y_{1}]\|=O\big{(}h^{p+\frac{1}{2}}\big{)}$
and, by Jensen’s inequality, we have

$\|\mathbb{E}[R_{p}(h,y_{0})]\|\leq\mathbb{E}\big{[}\|R_{p}(h,y_{0})\|^{2}\big{]}^{\frac{1}{2}},\hskip
8.53581pt\text{and}\hskip
8.53581pt\|\mathbb{E}[R_{p}^{\gamma}(h,y_{0})]\|\leq\mathbb{E}\big{[}\|R_{p}^{\gamma}(h,y_{0})\|^{2}\big{]}^{\frac{1}{2}}\hskip
0.7113pt.$

Since the above terms were previously shown to be $O(h^{p+\frac{1}{2}})$, the
result follows.

As a consequence of the local estimates given by Theorem 3.8 and following the
mean-square framework of Milstein and Tretyakov [59], we now obtain our main
result.

######  Theorem 3.10 (Global strong error estimate).

Given a fixed number of steps $N$, we define a numerical solution
$\\{Y_{k}\\}_{0\hskip 0.7113pt\leq\hskip 0.7113ptk\hskip 0.7113pt\leq\hskip
0.7113ptN}$ for the SDE (1) over $[0,T]$ as follows,

$Y_{k+1}:=\big{(}\text{Solution at }r=1\text{ of CDE (\ref{eq:intro_CDE})
driven by }\gamma_{k}:[0,1]\rightarrow\mathbb{R}^{1+d}\hskip
0.7113pt\big{)}\big{(}Y_{k}\big{)}+E_{k}\hskip 0.7113pt,$

where $Y_{0}:=y_{0}$ and, for a fixed
$p\in\\{\frac{k}{2}\\}_{k\in\mathbb{N}}\hskip 0.7113pt$, the “CDE errors”
$\\{E_{k}\\}$ uniformly satisfy

$\displaystyle\,\hskip
0.7113pt\mathbb{E}\big{[}\|E_{k}\|\big{]}^{\frac{1}{2}}=O(h^{p}),$
$\displaystyle\big{\|}\mathbb{E}[E_{k}]\big{\|}=O\big{(}h^{p+\frac{1}{2}}\big{)},$

and each path $\gamma_{k}:[0,1]\rightarrow\mathbb{R}^{1+d}$ is expressible as
$\gamma_{k}=\varphi\big{(}\big{\\{}(t,W_{t}):t\in\big{[}\frac{kT}{N},\frac{(k+1)T}{N}\big{]}\big{\\}}\big{)}$
for some fixed path-valued function $\varphi$. We will assume that the paths
$\\{\gamma_{k}\\}$ uniformly satisfy Assumption A1 and that
$f\in\mathcal{C}_{\mathrm{Lip}}^{\lceil
p-1\rceil}(\mathbb{R}^{e},\mathbb{R}^{e})$ and
$g\in\mathcal{C}_{\mathrm{Lip}}^{2p-1}(\mathbb{R}^{e},\mathbb{R}^{e\times
d})$. Suppose also that the assumptions of Theorem 2.2 hold and that the
integrals $I^{\gamma_{k}}_{\alpha}(1)$ and $I_{\alpha}(1)$ agree almost surely
for all $\alpha\in\mathcal{A}$ with $\textrm{ord}(\alpha)\leq p-\frac{1}{2}$
and in expectation for all $\alpha\in\mathcal{A}$ with
$\textrm{ord}(\alpha)=p$. Then over the finite interval $[0,T]$, for
$k\in\\{1,2,\cdots,N\\}$, we have

$\mathbb{E}\big{[}\hskip
0.7113pt\|y_{kh}-Y_{k}\|^{2}\big{]}^{1/2}=O\big{(}h^{p-\frac{1}{2}}\big{)}.$

######  Remark 3.11 (Infinite time horizon).

The convergence results in this paper are established over a finite time
horizon $T$. However, our framework could be employed to deal with the
infinite time horizon setting under suitable conditions on the SDE. For
instance, if the SDE is ergodic with an exponential contraction property then
contributions of local errors to the global error are reduced (exponentially
in time). An extension of the classical Milstein-Tretyakov mean-square error
analysis [59] to the infinite time horizon case for such contractive SDEs is
given by [52, Theorem 3.3.]. See also [25, Section 3] for similar, but
employing Multilevel Monte Carlo (MLMC).

###  3.4 Application to commutative SDEs

Although Theorem 3.10 identifies conditions on the piecewise linear path
$\gamma$ to achieve a given strong convergence rate, it can be difficult to
generate the required integrals (as discussed in the introduction).
Fortunately, the commutativity condition (2) leads to certain simplifications
in the Taylor expansions of the Stratonovich SDE (1) and its CDE approximation
(3). This is detailed in Appendix A but, as a consequence, we have the
following theorem.

######  Theorem 3.12 (Global strong error estimate for commutative SDEs).

Suppose that the assumptions of Theorem 3.10 hold for $p=2$, but with the
exception that each path $\gamma_{k}$ will now match only the following
iterated integrals of the Brownian motion:

$\displaystyle\gamma^{\tau}_{k}(1)-\gamma^{\tau}_{k}(0)$ $\displaystyle=h,$
$\displaystyle\gamma^{\hskip 0.7113pt\omega}_{k}(1)-\gamma^{\hskip
0.7113pt\omega}_{k}(0)$ $\displaystyle=W_{kh,(k+1)h}\hskip 0.7113pt,$
$\displaystyle\int_{0}^{1}\big{(}\gamma_{k}^{\hskip
0.7113pt\omega}(r)-\gamma_{k}^{\hskip 0.7113pt\omega}(0)\big{)}\hskip
0.7113ptd\gamma_{k}^{\tau}(r)$
$\displaystyle=\int_{kh}^{(k+1)h}W_{kh,u}\,du\hskip 0.7113pt,$
$\displaystyle\mathbb{E}\bigg{[}\int_{0}^{1}\Big{(}\big{(}\gamma_{k}^{\hskip
0.7113pt\omega}\big{)}^{i}(r)-\big{(}\gamma_{k}^{\hskip
0.7113pt\omega}\big{)}^{i}(0)\Big{)}^{2}d\gamma_{k}^{\tau}(r)\bigg{]}$
$\displaystyle=\mathbb{E}\bigg{[}\int_{kh}^{(k+1)h}W_{kh,u}^{2}\,du\hskip
0.7113pt\bigg{]}=\frac{1}{2}\hskip 0.7113pth^{2},$

and for $i,j\in\\{1,\cdots,d\\}$ with $i\neq j$, we have

$\displaystyle\mathbb{E}\bigg{[}\int_{0}^{1}\int_{0}^{r_{1}}\int_{0}^{r_{2}}d\big{(}\gamma_{k}^{\hskip
0.7113pt\omega}\big{)}^{i}(r_{3})\,d\big{(}\gamma_{k}^{\hskip
0.7113pt\omega}\big{)}^{j}(r_{2})\,d\gamma_{k}^{\tau}(r_{1})\bigg{]}=0,$
$\displaystyle\mathbb{E}\bigg{[}\int_{0}^{1}\int_{0}^{r_{1}}\int_{0}^{r_{2}}d\big{(}\gamma_{k}^{\hskip
0.7113pt\omega}\big{)}^{i}(r_{3})\,d\gamma_{k}^{\tau}(r_{2})\,d\big{(}\gamma_{k}^{\hskip
0.7113pt\omega}\big{)}^{j}(r_{1})\bigg{]}=0,$
$\displaystyle\mathbb{E}\bigg{[}\int_{0}^{1}\int_{0}^{r_{1}}\int_{0}^{r_{2}}d\gamma_{k}^{\tau}(r_{3})\,d\big{(}\gamma_{k}^{\hskip
0.7113pt\omega}\big{)}^{i}(r_{2})\,d\big{(}\gamma_{k}^{\hskip
0.7113pt\omega}\big{)}^{j}(r_{1})\bigg{]}=0.$

Then on the interval $[0,T]$, for $k\in\\{1,2,\cdots,N\\}$, the numerical
solution $\\{Y_{k}\\}$ satisfies

$\mathbb{E}\big{[}\,\|y_{kh}-Y_{k}\|^{2}\big{]}^{1/2}=O\big{(}h^{\frac{3}{2}}\big{)}.$

######  Proof 3.13.

By Theorem A.9, we see that the CDE Taylor expansion can match all the “noise
only” terms with $\textrm{ord}(\alpha)\leq 2$ simply by the path $\gamma_{k}$
having the increment $\gamma^{\hskip 0.7113pt\omega}_{k}(1)-\gamma^{\hskip
0.7113pt\omega}_{k}(0)=W_{kh,(k+1)h}\hskip 0.7113pt$. Adopting the general
notation used within Appendix A, we recall Theorem A.5, which gives the
following decompositions of iterated integrals:

(22) $\displaystyle I_{ij}$ $\displaystyle=\frac{1}{2}I_{i}\cdot
I_{j}+\frac{1}{2}I_{[i,j]}\hskip 0.7113pt,$ (23) $\displaystyle I_{ijk}$
$\displaystyle=\frac{1}{6}I_{i}\cdot I_{j}\cdot I_{k}+\frac{1}{4}I_{i}\cdot
I_{[j,k]}+\frac{1}{4}I_{[i,j]}\cdot
I_{k}+\frac{1}{6}\big{(}I_{[[i,j],k]}+I_{[i,[j,k]]}),$

for indices $i,j,k\in\\{1,\cdots,d\\}$. However, by identifying a coordinate
of the Brownian motion with time, we see that the above would still hold when
$i,j,k\in\\{0,1,\cdots,d\\}$. (In the following paragraph, the verb “match”
refers to when $r\mapsto(\gamma_{k}^{\tau},\gamma_{k}^{\hskip
0.7113pt\omega})(r)$ and $t\mapsto(t,W_{t})\hskip 0.7113pt$ give the same
iterated integral – either almost surely or in expectation).

Thus, by virtue of matching $\\{I_{0}\hskip 0.7113pt,I_{i}\hskip
0.7113pt,I_{i0}\\}_{1\hskip 0.7113pt\leq\hskip 0.7113pti\hskip
0.7113pt\leq\hskip 0.7113ptd}\hskip 0.7113pt$, $\gamma_{k}$ matches
$\\{I_{[i,0]}\\}$ and thus $\\{I_{0i}\\}\hskip 0.7113pt$. Using integration by
parts, we note the following identities for triple iterated integrals:

$\displaystyle\int_{0}^{1}\int_{0}^{r_{1}}\int_{0}^{r_{2}}d\big{(}\gamma_{k}^{\hskip
0.7113pt\omega}\big{)}^{i}(r_{3})\,d\big{(}\gamma_{k}^{\hskip
0.7113pt\omega}\big{)}^{i}(r_{2})\,d\gamma_{k}^{\tau}(r_{1})$
$\displaystyle=\int_{0}^{1}\frac{1}{2}\Big{(}\big{(}\gamma_{k}^{\hskip
0.7113pt\omega}\big{)}^{i}(r)-\big{(}\gamma_{k}^{\hskip
0.7113pt\omega}\big{)}^{i}(0)\Big{)}^{2}d\gamma_{k}^{\tau}(r),$
$\displaystyle\int_{kh}^{(k+1)h}\int_{kh}^{u_{1}}\int_{kh}^{u_{2}}\circ\,dW_{u_{3}}^{i}\hskip
0.7113pt\circ\hskip 0.7113ptdW_{u_{2}}^{i}\,du_{1}$
$\displaystyle=\int_{kh}^{(k+1)h}\frac{1}{2}\hskip
0.7113ptW_{kh,u}^{2}\,du\hskip 0.7113pt.$

Since $\gamma_{k}$ is assumed to match $\\{I_{ii0}\\}$ in expectation and the
lower order terms exactly, by (23) and the fact that $I_{[0,[i,i]]}=0$, it
will also match $\\{I_{[i,[i,0]]}\\}$ in expectation. Hence, $\gamma_{k}$
matches $\\{I_{i0i},I_{0ii}\\}$ in expectation by (23) and the antisymmetry of
$[\,\cdot\hskip 0.7113pt,\hskip 0.7113pt\cdot\hskip 0.7113pt]$. By assumption,
the remaining integrals $\\{I_{ij0},I_{i0j},I_{0ij}\\}$ are matched in
expectation.

From the above, we see the Taylor expansions of the SDE (1) and CDE (3)
coincide up to order $p=2$, as required by Theorem 3.10. The result now
follows.

######  Remark 3.14.

In Theorem 3.10, we have accounted for the fact that the CDE (3), or rather
the resulting sequence of ODEs, may be approximated using an ODE solver.
However, obtaining the required estimates for these additional “CDE errors”
$\\{E_{k}\\}$ may be non-trivial and thus, we leave such an error analysis as
a topic of future work. That said, to achieve strong order $3/2$ convergence,
we expect that a single step of a second order ODE solver would be sufficient
to discretise ODEs depending on just $f$ and a single step of a fourth order
solver (such as RK4) to suffice for the other ODEs. The intuition is that
$\gamma$ has Brownian-like scaling and so vector fields are either $O(h)$ or
$O\big{(}h^{\frac{1}{2}}\big{)}$. Hence, we expect the local errors to be
$O(h^{3})$ or $O\big{(}h^{\frac{5}{2}}\big{)}$ in these two cases.

##  4 Paths

In this section, we present a variety of piecewise linear paths which fall
into the proposed framework for developing SDE splitting methods (Theorem
3.10). These “splitting paths” correspond to both well-known numerical methods
(such as Lie-Trotter and Strang splitting [77]) as well as the new high order
splitting methods, which can exploit the optimal integral estimators that are
derived in Appendix B. Furthermore, we illustrate both the Strang and high
order splitting paths in Figure 3. Throughout, we use the notation in Example
1.3 and define paths by their increments.

######  Example 4.1 (Lie-Trotter).

A Lie-Trotter splitting can be defined by one of two possible two-piece paths
$\gamma^{LT1},\gamma^{LT2}:[0,1]\rightarrow\mathbb{R}^{1+d}$ given by
$\gamma^{LT}(z)=(\gamma^{\tau},\gamma^{\hskip 0.7113pt\omega})(z)$ with

(24) $\displaystyle\gamma_{r_{i},r_{i+1}}^{LT1}$
$\displaystyle:=\begin{cases}(h,0),\quad&\text{if }\,i=0\\\\[6.0pt]
(0,W_{s,t}),&\text{if }\,i=1,\end{cases}$ (25)
$\displaystyle\gamma_{r_{i},r_{i+1}}^{LT2}$
$\displaystyle:=\begin{cases}\left(0,W_{s,t}\right),\quad&\text{if
}\,i=0\\\\[6.0pt] \left(h,0\right),&\text{if }\,i=1.\end{cases}$

######  Example 4.2 (Strang splitting).

The Strang splitting, see Figure 3, can be defined as a three-piece path
$\gamma^{S}:[0,1]\rightarrow\mathbb{R}^{1+d}$ given by
$\gamma^{S}(z)=(\gamma^{\tau},\gamma^{\hskip 0.7113pt\omega})(z)$ with the
pieces:

(26)
$\displaystyle\gamma_{r_{i},r_{i+1}}^{S}:=\begin{cases}\big{(}\frac{1}{2}h,0\big{)},&\text{if
}\,i=0\\\\[6.0pt] (0,W_{s,t}),&\text{if }\,i=1\\\\[6.0pt]
\big{(}\frac{1}{2}h,0\big{)},&\text{if }\,i=2.\end{cases}$

We now proceed to “higher order” piecewise linear paths that are constructed
to match the increment $W_{s,t}$ and space-time Lévy area $H_{s,t}$ of the
Brownian motion. For these paths to match higher order iterated integrals in
expectation (for example, to achieve 3/2 strong convergence), then we
necessarily require at least three pieces. The inability of paths with two
pieces to match the conditions (8) and (9) required for high order strong
convergence was explicitly shown in [26, p97 and Appendix A]. We begin by
presenting paths (27) and (28), which each have a total of five pieces
(vertical and horizontal), and can thus be seen as extensions of the Strang
splitting.

######  Example 4.3 (High order Strang splitting (linear version)).

A high order Strang splitting, see Figure 3, can be defined using a five-piece
path $\gamma^{HS1}:[0,1]\rightarrow\mathbb{R}^{1+d}$, which is linear in the
Brownian motion and has the pieces:

(27)
$\displaystyle\gamma_{r_{i},r_{i+1}}^{HS1}:=\begin{cases}\big{(}\frac{3-\sqrt{3}}{6}h,0\big{)},&\text{if
}\,i=0\\\\[6.0pt] \big{(}0,\frac{1}{2}W_{s,t}+\sqrt{3}H_{s,t}\big{)},&\text{if
}\,i=1\\\\[6.0pt] \big{(}\frac{\sqrt{3}}{3}h,0\big{)},&\text{if
}\,i=2\\\\[6.0pt] \big{(}0,\frac{1}{2}W_{s,t}-\sqrt{3}H_{s,t}\big{)},&\text{if
}\,i=3\\\\[6.0pt] \big{(}\frac{3-\sqrt{3}}{6}h,0\big{)},&\text{if
}\,i=4.\end{cases}$

######  Example 4.4 (High order Strang splitting (non-linear version)).

A high order Strang splitting, see Figure 3, can be defined as a five-piece
path $\gamma^{HS2}:[0,1]\rightarrow\mathbb{R}^{1+d}$, which is based on an
optimal estimator for a certain Brownian integral and has pieces:

(28)
$\displaystyle\gamma_{r_{i},r_{i+1}}^{HS2}:=\begin{cases}\big{(}0,\frac{1}{2}W_{s,t}+H_{s,t}-\frac{1}{2}C_{s,t}\big{)},&\text{if
}\,i=0\\\\[5.0pt] \big{(}\frac{1}{2}h,0\big{)},&\text{if }\,i=1\\\\[5.0pt]
\big{(}0,C_{s,t}\big{)},&\text{if }\,i=2\\\\[5.0pt]
\big{(}\frac{1}{2}h,0\big{)},&\text{if }\,i=3\\\\[5.0pt]
\big{(}0,\frac{1}{2}W_{s,t}-H_{s,t}-\frac{1}{2}C_{s,t}\big{)},&\text{if
}\,i=4.\end{cases}$

where the random vector $C_{s,t}$ is defined component-wise by

(29) $\displaystyle C_{s,t}^{j}$
$\displaystyle:=\epsilon_{s,t}^{j}\bigg{(}\frac{1}{3}\big{(}W_{s,t}^{j}\big{)}^{2}+\frac{4}{5}\big{(}H_{s,t}^{j}\big{)}^{2}+\frac{4}{15}h-\frac{1}{\sqrt{6\pi}}h^{\frac{1}{2}}n_{s,t}^{j}W_{s,t}^{j}\bigg{)}^{\frac{1}{2}},$
(30) $\displaystyle\epsilon_{s,t}^{j}$
$\displaystyle:=\mathrm{sgn}\bigg{(}W_{s,t}^{j}-\frac{3}{\sqrt{24\pi}}h^{\frac{1}{2}}n_{s,t}^{j}\bigg{)},$

where
$n_{s,t}^{j}:=\mathrm{sgn}\big{(}H_{s,s+\frac{1}{2}h}^{j}-H_{s+\frac{1}{2}h,t}^{j}\big{)}$
are independent Rademacher random variables.

######  Remark 4.5.

The formula for $C_{s,t}$ is derived such that
$\int_{0}^{1}\big{(}\big{(}\gamma_{r}^{\hskip
0.7113pt\omega}\big{)}^{j}-\big{(}\gamma_{0}^{\hskip
0.7113pt\omega}\big{)}^{j}\hskip 0.7113pt\big{)}^{2}d\gamma_{r}^{\tau}$
matches the optimal estimator
$\mathbb{E}\big{[}\int_{s}^{t}\big{(}W_{s,u}^{j}\big{)}^{2}\,du\,\big{|}\,W_{s,t}\hskip
0.7113pt,H_{s,t}\hskip 0.7113pt,n_{s,t}\big{]}$ (see Theorem B.9).

For paths with three pieces, two of which are vertical and only relate to
diffusion vector field, we refer to the resulting approximation as the
“Shifted ODE” approach. We use this terminology as, in the additive noise
setting, the vertical pieces correspond to additive shifts for the numerical
solution and so there is only one non-trivial ODE. As before, paths can be
linear or non-linear functions of the input random variables.

######  Example 4.6 (Shifted ODE splitting (high order and linear)).

We can define a high order splitting by a three-piece path
$\gamma^{SO1}:[0,1]\rightarrow\mathbb{R}^{1+d}$ with the following pieces:

(31)
$\displaystyle\gamma_{r_{i},r_{i+1}}^{SO1}:=\begin{cases}\big{(}0,H_{s,t}+\frac{1}{2}\sqrt{h}n_{s,t}\big{)},&\text{if
}\,i=0\\\\[5.0pt] \big{(}h,W_{s,t}-\sqrt{h}n_{s,t}\big{)},&\text{if
}\,i=1\\\\[5.0pt]
\big{(}0,-H_{s,t}+\frac{1}{2}\sqrt{h}n_{s,t}\big{)},&\text{if
}\,i=2.\end{cases}$

######  Example 4.7 (Shifted ODE splitting (high order and non-linear)).

We can define a high order splitting, see Figure 3, using a three-piece path
$\hskip 0.7113pt\gamma^{SO2}:[0,1]\rightarrow\mathbb{R}^{1+d}$, which is based
on an optimal estimator for a certain Brownian integral and has pieces:

(32)
$\displaystyle\gamma_{r_{i},r_{i+1}}^{SO2}:=\begin{cases}\big{(}0,\frac{1}{2}W_{s,t}+H_{s,t}-\frac{1}{2}C_{s,t}\big{)},&\text{if
}\,i=0\\\\[5.0pt] \big{(}h,C_{s,t}\big{)},&\text{if }\,i=1\\\\[5.0pt]
\big{(}0,\frac{1}{2}W_{s,t}-H_{s,t}-\frac{1}{2}C_{s,t}\big{)},&\text{if
}\,i=2.\end{cases}$

where the random vector $C_{s,t}$ is defined component-wise by

$\displaystyle C_{s,t}^{j}$
$\displaystyle:=\epsilon_{s,t}^{j}\bigg{(}\big{(}W_{s,t}^{j}\big{)}^{2}+\frac{12}{5}\big{(}H_{s,t}^{j}\big{)}^{2}+\frac{4}{5}h-\frac{3}{\sqrt{6\pi}}h^{\frac{1}{2}}n_{s,t}^{j}W_{s,t}^{j}\bigg{)}^{\frac{1}{2}},$
$\displaystyle\epsilon_{s,t}^{j}$
$\displaystyle:=\mathrm{sgn}\bigg{(}W_{s,t}^{j}-\frac{3}{\sqrt{24\pi}}h^{\frac{1}{2}}n_{s,t}^{j}\bigg{)},$

where
$n_{s,t}^{j}:=\mathrm{sgn}\big{(}H_{s,s+\frac{1}{2}h}^{j}-H_{s+\frac{1}{2}h,t}^{j}\big{)}$
are independent Rademacher random variables.

######  Remark 4.8.

Just like Example 4.4, $C_{s,t}$ is derived so that
$\int_{0}^{1}\big{(}\big{(}\gamma_{r}^{\hskip
0.7113pt\omega}\big{)}^{j}-\big{(}\gamma_{0}^{\hskip
0.7113pt\omega}\big{)}^{j}\hskip 0.7113pt\big{)}^{2}d\gamma_{r}^{\tau}$
matches the optimal estimator
$\mathbb{E}\big{[}\int_{s}^{t}\big{(}W_{s,u}^{j}\big{)}^{2}\,du\,\big{|}\,W_{s,t}\hskip
0.7113pt,H_{s,t}\hskip 0.7113pt,n_{s,t}\big{]}$ (see Theorem B.9).

The following paths do not generally result in high order approximations for
SDEs satisfying the commutativity condition (2). However, the piecewise linear
path given by (33) results in the “Shifted Euler” method for SDEs with
additive noise, which we demonstrate can outperform the standard Euler-
Maruyama method in Section 5.

######  Example 4.9 (Shifted ODE splitting (low order; suitable for Euler’s
method)).

We can define a low order splitting by a three-piece path
$\gamma^{SO3}:[0,1]\rightarrow\mathbb{R}^{1+d}$ with

(33)
$\displaystyle\gamma_{r_{i},r_{i+1}}^{SO3}:=\begin{cases}\big{(}0,\frac{1}{2}W_{s,t}+H_{s,t}\big{)},&\text{if
}\,i=0\\\\[6.0pt] \big{(}h,0\big{)},&\text{if }\,i=1\\\\[6.0pt]
\big{(}0,\frac{1}{2}W_{s,t}-H_{s,t}\big{)},&\text{if }\,i=2.\end{cases}$

The path (34) is also not usually high order, but gives a third order
approximation when applied to underdamped Langevin dynamics (see Section 5 or
[30] for details).

######  Example 4.10 (Shifted ODE splitting for the underdamped Langevin
diffusion [30]).

We can define a splitting by a three-piece path
$\gamma^{SO4}:[0,1]\rightarrow\mathbb{R}^{1+d}$ with pieces:

(34)
$\displaystyle\gamma_{r_{i},r_{i+1}}^{SO4}:=\begin{cases}\big{(}0,H_{s,t}+6K_{s,t}\big{)},&\text{if
}\,i=0\\\\[6.0pt] \big{(}h,W_{s,t}-12K_{s,t}\big{)},&\text{if
}\,i=1\\\\[6.0pt] \big{(}0,-H_{s,t}+6K_{s,t}\big{)},&\text{if
}\,i=2.\end{cases}$

where $K_{s,t}\sim\mathcal{N}\big{(}0,\frac{1}{720}hI_{d}\big{)}$ is
independent of $\big{(}W_{s,t}\hskip 0.7113pt,H_{s,t}\big{)}$ and given by the
integral

$\displaystyle K_{s,t}$
$\displaystyle:=\frac{1}{h^{2}}\int_{s}^{t}\bigg{(}W_{s,u}-\frac{u-s}{h}\,W_{s,t}\bigg{)}\bigg{(}\frac{1}{2}h-(u-s)\bigg{)}du\hskip
0.7113pt.$

Figure 3: Illustration of piecewise linear paths associated with various
splitting methods for SDEs.

##  5 Examples

We demonstrate the proposed splitting methods on several SDEs (the Cox-
Ingersoll-Ross model and additive noise SDEs such as a stochastic oscillator,
the stochastic FitzHugh-Nagumo model and the underdamped Langevin diffusion).
For the majority of these examples, the “non-diffusion” ODEs coming from the
SDE splitting will not admit a closed-formed solution and thus must be further
discretised. We will show that such ODEs can be resolved by Runge-Kutta or
splitting methods.

Throughout, we shall compare methods using the following strong error
estimator:

######  Definition 5.1 (Strong error estimator for SDEs).

For $N\geq 1$, let $Y_{N}$ denote a numerical solution to the SDE (1) computed
at time $T$ with a fixed step size $h=\frac{T}{N}\hskip 0.7113pt$. Then we
define the following estimator for quantifying the strong convergence of
$Y_{N}$:

(35)
$S_{N}:=\sqrt{\mathbb{E}\Big{[}\big{(}Y_{N}-Y_{T}^{\text{fine}}\big{)}^{2}\Big{]}},$

where $Y_{T}^{\text{fine}}$ denotes a numerical solution to (1) computed with
a finer step size, $h^{\text{fine}}\leq\frac{1}{10}\hskip 0.7113pth$, but
using the same Brownian motion (so that $Y_{N}$ and $Y_{T}^{\text{fine}}$ are
close). In our examples, the expectation in (35) will be estimated by standard
Monte Carlo.

All the experiments were conducted on a laptop in either C++, Python or R.
Associated code is available at github.com/james-m-foster/high-order-splitting
and github.com/james-m-foster/high-order-langevin for the Langevin dynamics
example.

###  5.1 Cox-Ingersoll-Ross model

The Cox-Ingersoll-Ross (or CIR) model is a popular one-factor short rate model
used in mathematical finance for modelling interest rates [17] and stochastic
volatility [38]. It is given by the following Itô SDE:

(36) $dy_{t}=a(b-y_{t})\,dt+\sigma\sqrt{y_{t}}\,dW_{t}\hskip 0.7113pt,$

where the parameters $a,b,\sigma\geq 0$ describe the mean reversion
speed/level and volatility. Over the years, a variety of numerical methods
have been proposed for the CIR model and we refer the reader to some of these
approaches [2, 3, 4, 18, 22, 26, 37, 42, 58, 67]. As we propose splitting
methods, we first rewrite the SDE (36) in Stratonovich form:

$dy_{t}=a(\hskip 0.7113pt\widetilde{b}-y_{t})\,dt+\sigma\sqrt{y_{t}}\circ
dW_{t}\hskip 0.7113pt.$

where $\widetilde{b}:=b-\frac{\sigma^{2}}{4a}\hskip 0.7113pt$. To ensure non-
negativity of $\widetilde{b}$ and our scheme, we assume $\sigma^{2}\leq 4ab$.
By driving (36) with the piecewise linear path (27), we obtain the splitting
method:

$\displaystyle Y_{k}^{(1)}$
$\displaystyle:=e^{-\frac{3-\sqrt{3}}{6}ah}Y_{k}+\widetilde{b}\big{(}1-e^{-\frac{3-\sqrt{3}}{6}ah}\big{)},$
$\displaystyle Y_{k}^{(2)}$
$\displaystyle:=\bigg{(}\sqrt{Y_{k}^{(1)}}+\frac{\sigma}{2}\Big{(}\frac{1}{2}W_{k}+\sqrt{3}H_{k}\Big{)}\bigg{)}^{2},$
$\displaystyle Y_{k}^{(3)}$
$\displaystyle:=e^{-\frac{\sqrt{3}}{3}ah}Y_{k}^{(2)}+\widetilde{b}\big{(}1-e^{-\frac{\sqrt{3}}{3}ah}\big{)},$
$\displaystyle Y_{k}^{(4)}$
$\displaystyle:=\bigg{(}\sqrt{Y_{k}^{(3)}}+\frac{\sigma}{2}\Big{(}\frac{1}{2}W_{k}-\sqrt{3}H_{k}\Big{)}\bigg{)}^{2},$
(37) $\displaystyle Y_{k+1}$
$\displaystyle:=e^{-\frac{3-\sqrt{3}}{6}ah}Y_{k}^{(4)}+\widetilde{b}\big{(}1-e^{-\frac{3-\sqrt{3}}{6}ah}\big{)},$

since for the CIR model, the drift and diffusion ODEs admit closed-form
solutions. The resulting numerical method (5.1) is then straightforward to
implement and we might expect each step to be roughly twice as expensive as
previous methodologies (since the above method requires generating two
Gaussian random variables per step). We also note that the method (5.1), and
each of its stages, preserves non-negativity. In addition, one can compute the
mean and variance of both the SDE (36) and (5.1).

######  Theorem 5.2.

The numerical solution given by (5.1) has the following moments:

$\displaystyle\mathbb{E}[Y_{k+1}|Y_{k}]$
$\displaystyle=e^{-ah}Y_{k}+b\big{(}1-e^{-ah}\big{)}+R_{k}^{E}\hskip
0.7113pt,$ $\displaystyle\operatorname{Var}(Y_{k+1}|Y_{k})$
$\displaystyle=\frac{\sigma^{2}}{a}\big{(}e^{-ah}-e^{-2ah}\big{)}Y_{k}+\frac{b\sigma^{2}}{2a}\big{(}1-e^{-ah}\big{)}^{2}+R_{k}^{V}\hskip
0.7113pt,$

where the remainder terms $R_{k}^{E}$ and $R_{k}^{V}$ can be estimated as
$O(h^{5})$ in an $L^{2}(\mathbb{P})$ sense.

######  Proof 5.3.

See the proof of Theorem C.1 in the appendix.

######  Remark 5.4.

The associated mean and variance of the CIR model (36) are also known
analytically [17] and are given by the above formulae without the $O(h^{5})$
terms.

Therefore, even though the primary focus of the paper is strong approximation,
we see that the splitting (5.1) is also a high order weak approximation of the
CIR model. Furthermore, due to the non-Lipschitz diffusion vector field in the
CIR model, we can not directly apply our main result to establish the strong
convergence of our method. We thus leave the error analysis of path-based
splitting methods in the “non-smooth” setting as future work. Similar to [8],
we expect results can be obtained in cases where the drift has polynomial
growth and satisfies a (global) one-sided Lipschitz condition.

We will now present our experiment, comparing the strong convergence of the
proposed splitting (5.1) against standard methods. We use the following
parameters:

$a=1,\hskip 14.22636ptb=1,\hskip 14.22636pt\sigma=1,\hskip
14.22636pty_{0}=1,\hskip 14.22636ptT=1.$

Figure 4: $S_{N}$ estimated for (36) with 100,000 sample paths as a function
of step size $h=\frac{T}{N}$.

In the above graph, we see that the proposed splitting scheme achieves a high
order convergence rate and better accuracy than other schemes (for a fixed
step size). Moreover, we believe (5.1) is the first non-negative method for
the CIR model to show high order strong convergence and to match the first two
moments with $O(h^{5})$ error. However, since the square root vector field is
not globally Lipschitz continuous, our scheme converges slightly slower than
the $O(h^{\frac{3}{2}})$ rate obtained in the “smooth” case.

As we discussed previously, each step of the splitting method (5.1) will be
more computationally expensive than for low order schemes. This is quantified
in Table 1.

Table 1: Computer time to simulate $100,000$ sample paths of (36) in C++ with $100$ steps (seconds) Splitting (5.1) | Ninomiya-Victoir | Drift-Implicit Euler | Milstein | Euler
---|---|---|---|---
2.13 | 1.07 | 1.42 | 1.01 | 0.86

We see that the high order splitting method (5.1) is roughly twice as
expensive as the lower order schemes. Combining these times with Figure 4, we
arrive at Table 2, which shows our method achieves a small error significantly
faster than other schemes.

Table 2: Estimated time to produce $100,000$ sample paths of (36) with an error of $S_{N}=10^{-3}$ (seconds) Splitting (5.1) | Ninomiya-Victoir | Drift-Implicit Euler | Milstein | Euler
---|---|---|---|---
0.27 | 1.99 | 4.17 | 3.69 | 490

###  5.2 Additive noise SDEs, such as a stochastic anharmonic oscillator

We first consider an arbitrary SDE with additive noise (that is, $g(\cdot)$ is
a fixed matrix)

(38) $dy_{t}=f(y_{t})\,dt+\sigma\hskip 0.7113ptdW_{t}\hskip 0.7113pt,$

where $\sigma\in\mathbb{R}^{e\times d}$. As (38) has additive noise, it is in
both Itô and Stratonovich form. Since each piece in $\gamma$ yields an ODE, we
would like to minimize the number of pieces. Moreover, it was shown in [26,
p97 and Appendix A] that two -piece paths are unable to match the various
iterated integrals of Brownian motion required in Theorem 3.12. Thus, in
general, we propose driving (38) by a piecewise linear path with three pieces.
So unless the SDE has special structure, such as underdamped Langevin
dynamics, we recommend either the paths (31) or (32) to achieve 3/2 order
strong convergence. In either case, we obtain a splitting method for the SDE
(38) with the following form:

(39) $Y_{k+1}:=\exp\big{(}f(\cdot)h+\sigma\hskip
0.7113ptC_{2}\big{)}\big{(}Y_{k}+\sigma\hskip
0.7113ptC_{1}\big{)}+\sigma\hskip 0.7113ptC_{3}\hskip 0.7113pt,$

where the vectors $C_{1}\hskip 0.7113pt,C_{2}\hskip
0.7113pt,C_{3}\in\mathbb{R}^{d}$ correspond to the increments of the driving
path. Since the ODE in (39) may not admit a closed-form solution, it must be
discretised. Expanding terms in the Taylor expansion of this ODE, we see that
(39) is equal to

$\displaystyle\bigg{(}I_{d}(\cdot)+f(\cdot)h+\sigma C_{2}+\frac{1}{2}f^{\hskip
0.7113pt\prime}(\cdot)\big{(}f(\cdot)h+\sigma
C_{2}\big{)}h+\frac{1}{6}f^{\hskip 0.7113pt\prime\prime}(\cdot)\big{(}\sigma
C_{2}\big{)}^{\otimes 2}h\bigg{)}\big{(}Y_{k}+\sigma\hskip
0.7113ptC_{1}\big{)}+\sigma\hskip 0.7113ptC_{3}$
$\displaystyle=Y_{k}+\sigma\hskip
0.7113ptC_{1}+\bigg{(}f\big{(}Y_{k}\big{)}+f^{\hskip
0.7113pt\prime}\big{(}Y_{k}\big{)}\sigma\hskip
0.7113ptC_{1}+\frac{1}{2}f^{\hskip
0.7113pt\prime\prime}\big{(}Y_{k}\big{)}\big{(}\sigma\hskip
0.7113ptC_{1}\big{)}^{\otimes 2}\bigg{)}h+\sigma\hskip
0.7113ptC_{2}+\frac{1}{2}f^{\hskip
0.7113pt\prime}\big{(}Y_{k}\big{)}f\big{(}Y_{k}\big{)}h^{2}$
$\displaystyle\hskip 14.22636pt+\frac{1}{2}\Big{(}f^{\hskip
0.7113pt\prime}\big{(}Y_{k}\big{)}\sigma\hskip 0.7113ptC_{2}+f^{\hskip
0.7113pt\prime\prime}\big{(}Y_{k}\big{)}\big{(}\big{(}\sigma\hskip
0.7113ptC_{2}\big{)}\otimes\big{(}\sigma\hskip
0.7113ptC_{1}\big{)}\big{)}\Big{)}h+\frac{1}{6}f^{\hskip
0.7113pt\prime\prime}\big{(}Y_{k}\big{)}\big{(}\sigma C_{2}\big{)}^{\otimes
2}h+\sigma\hskip 0.7113ptC_{3}+\cdots$

Rearranging these terms then gives the following expansion for the splitting
method,

$\displaystyle Y_{k+1}\approx Y_{k}$
$\displaystyle+f\big{(}Y_{k}\big{)}h+\sigma\big{(}C_{1}+C_{2}+C_{3}\big{)}+\frac{1}{2}f^{\hskip
0.7113pt\prime}\big{(}Y_{k}\big{)}\bigg{(}\sigma\hskip
0.7113ptC_{1}+\frac{1}{2}\sigma\hskip 0.7113ptC_{2}\bigg{)}h$
$\displaystyle+\frac{1}{2}\hskip 0.7113ptf^{\hskip
0.7113pt\prime}\big{(}Y_{k}\big{)}f\big{(}Y_{k}\big{)}h^{2}+\frac{1}{2}\hskip
0.7113ptf^{\hskip
0.7113pt\prime\prime}\big{(}Y_{k}\big{)}\bigg{(}\big{(}\sigma\hskip
0.7113ptC_{1}\big{)}^{2}+\big{(}\sigma\hskip
0.7113ptC_{1}\big{)}\otimes\big{(}\sigma\hskip
0.7113ptC_{2}\big{)}+\frac{1}{3}\big{(}\sigma\hskip
0.7113ptC_{2}\big{)}^{2}\bigg{)}h,$

where the final line follows as the second derivative $f^{\hskip
0.7113pt\prime\prime}(Y_{k})$ is symmetric and bilinear. Since it will become
necessary to approximate the ODE map
$x\mapsto\exp\big{(}f(\cdot)h+\sigma\hskip 0.7113ptC_{2}\big{)}x$, we would
like the Taylor expansion of the ODE solver to coincide with the above. For
example, if we apply an explicit two-stage second order Runge-Kutta method
(determined by a parameter $\alpha$), then this will result in the following
Taylor expansion:

$\displaystyle x^{\text{RK}}$
$\displaystyle=x+\Big{(}1-\frac{1}{2\alpha}\Big{)}\big{(}f(x)h+\sigma
C_{2}\big{)}+\frac{1}{2\alpha}\Big{(}f\big{(}x+\alpha\big{(}f(x)h+\sigma
C_{2}\big{)}\big{)}h+\sigma C_{2}\Big{)}$ $\displaystyle=x+f(x)h+\sigma
C_{2}+\frac{1}{2}f^{\hskip 0.7113pt\prime}(x)\big{(}f(x)h+\sigma
C_{2}\big{)}h+\frac{1}{4}\alpha\hskip 0.7113ptf^{\hskip
0.7113pt\prime\prime}(x)(\sigma\hskip 0.7113ptC_{2})^{\otimes 2}h+\cdots.$

Thus, the only explicit two-stage Runge-Kutta method matching the Taylor
expansion of the splitting method (39) is Ralston’s method [70] (which
corresponds to $\alpha=\frac{2}{3}$). Hence, we propose the following
numerical method for the additive noise SDE (38),

$\displaystyle\widetilde{Y}_{k}^{\text{SR}}$
$\displaystyle:=Y_{k}^{\text{SR}}+\sigma\hskip 0.7113ptC_{1}\hskip 0.7113pt,$
$\displaystyle\widetilde{Y}_{k+\frac{2}{3}}^{\text{SR}}$
$\displaystyle:=\widetilde{Y}_{k}^{\text{SR}}+\frac{2}{3}\Big{(}f\big{(}\widetilde{Y}_{k}^{\text{SR}}\big{)}h+\sigma\hskip
0.7113ptC_{2}\Big{)}\hskip 0.7113pt,$ (40) $\displaystyle Y_{k+1}^{\text{SR}}$
$\displaystyle:=Y_{k}^{\text{SR}}+\frac{1}{4}f\big{(}\widetilde{Y}_{k}^{\text{SR}}\big{)}h+\frac{3}{4}f\big{(}\widetilde{Y}_{k+\frac{2}{3}}^{\text{SR}}\big{)}h+\sigma\hskip
0.7113ptW_{k}\hskip 0.7113pt,$

where $C_{1}\hskip 0.7113pt,C_{2}\in\mathbb{R}^{d}$ are the first two
increments of the driving piecewise linear path $\gamma$. Although we shall
not formally present an error analysis, it is clear from the above Taylor
expansions that the method (5.2) achieves a 3/2 strong order convergence rate.
Moreover, we can extend this method to SDEs with a time-varying diffusion
$\sigma\equiv\sigma(t)$. For example, the terms $\sigma\hskip
0.7113ptC_{1}\hskip 0.7113pt,\sigma\hskip 0.7113ptC_{2}$ and $\sigma\hskip
0.7113ptW_{k}$ may be replaced by $\sigma(t_{k})\hskip 0.7113ptC_{1}$,
$\sigma(t_{k})\hskip 0.7113ptC_{2}$ and

(41)
$\int_{t_{k}}^{t_{k+1}}\sigma(r)\,dW_{r}=\frac{1}{2}\big{(}\sigma(t_{k})+\sigma(t_{k+1})\big{)}W_{k}+\big{(}\sigma(t_{k+1})-\sigma(t_{k})\big{)}H_{k}+O\big{(}h^{2.5}\big{)}\hskip
0.7113pt.$

Similarly, we consider the case where the splitting path (33) is applied to
the additive noise SDE (38) and the resulting ODE is discretised using Euler’s
method. This gives

(42)
$Y_{k+1}^{\text{SE}}:=Y_{k}^{\text{SE}}+f\Big{(}Y_{k}^{\text{SE}}+\sigma\Big{(}\frac{1}{2}W_{k}+H_{k}\Big{)}\Big{)}h+\sigma\hskip
0.7113ptW_{k}\hskip 0.7113pt.$

We refer to the numerical method (42) as the “Shifted Euler” discretisation of
(38). Whilst the Shifted Euler method will have the same first order of
convergence as the Euler-Maruyama method, we expect it to be more accurate
(for sufficiently small $h$). This is due to the shifted Euler method (42)
having the following Taylor expansion:

$Y_{k+1}^{\text{SE}}=Y_{k}^{\text{SE}}+f\big{(}Y_{k}^{\text{SE}}\big{)}h+\sigma\hskip
0.7113ptW_{k}+f^{\hskip
0.7113pt\prime}\big{(}Y_{k}^{\text{SE}}\big{)}\bigg{(}\sigma\underbrace{\bigg{(}\frac{1}{2}W_{k}+H_{k}\bigg{)}h}_{=\,\int_{t_{k}}^{t_{k+1}}W_{t_{k}\hskip
0.7113pt,r}\,dr}\,\bigg{)}+O(h^{2}).$

Thus, the Shifted Euler scheme has a local error of $O(h^{2})$ whereas the
Euler-Maruyama method produces a local error of $O(h^{\frac{3}{2}})$, which is
due to the integral highlighted above. Similar to (5.2), the shifted Euler
method can be extended to SDEs with time-varying diffusion coefficients by
replacing terms $\sigma\hskip 0.7113ptW_{k}$ and
$\sigma\int_{t_{k}}^{t_{k+1}}W_{t_{k}\hskip 0.7113pt,r}\,dr$ with
$\int_{t_{k}}^{t_{k+1}}\sigma(r)\hskip 0.7113ptdW_{r}$ and
$\int_{t_{k}}^{t_{k+1}}\int_{t_{k}}^{s}\sigma(r)\,dW_{r}\,ds=\sigma(t_{k})\int_{t_{k}}^{t_{k+1}}W_{t_{k}\hskip
0.7113pt,r}\,dr+\frac{1}{6}\big{(}\sigma(t_{k+1})-\sigma(t_{k})\big{)}W_{k}\hskip
0.7113pth^{2}+O(h^{\frac{7}{2}})$.

Moreover, in applications where evaluating $f$ is significantly faster than
generating Gaussian random vectors, we can remove the space-time Lévy area
$H_{k}$ from (42). The resulting “increment-only” Shifted Euler method will
then have the same cost per step as the traditional Euler-Maruyama method.
However, for sufficiently small $h$, it will still be more accurate as it
gives $\frac{1}{2}hW_{k}\approx\int_{t_{k}}^{t_{k+1}}W_{t_{k}\hskip
0.7113pt,r}\,dr$ in its Taylor expansion.

We now present our numerical example; a scalar stochastic anharmonic
oscillator.

(43) $\hskip 113.81102ptdy_{t}=\sin(y_{t})\,dt+dW_{t}\hskip 0.7113pt,\hskip
51.21495pt(y_{0}=1,\hskip 7.11317ptT=1)$

We compare our approaches (high order method (5.2) and low order method (42))
against the 3/2 strong order SRA1 scheme in [71] and the Euler-Maruyama
method:

$\displaystyle Y_{k+1}^{\text{A1}}$
$\displaystyle:=Y_{k}^{\text{A1}}+\frac{1}{3}\hskip
0.7113ptf\big{(}Y_{k}^{\text{A1}}\big{)}h+\frac{2}{3}\hskip
0.7113ptf\bigg{(}Y_{k}^{\text{A1}}+\frac{3}{4}\Big{(}f\big{(}Y_{k}^{\text{A1}}\big{)}h+\sigma(W_{k}+2H_{k})\Big{)}\bigg{)}h+\sigma\hskip
0.7113ptW_{k}\hskip 0.7113pt,$ $\displaystyle Y_{k+1}^{\text{EM}}$
$\displaystyle:=Y_{k}^{\text{EM}}+f\big{(}Y_{k}^{\text{EM}}\big{)}h+\sigma\hskip
0.7113ptW_{k}\hskip 0.7113pt.$

We use the non-linear splitting path (32) to obtain the variables $C_{1}\hskip
0.7113pt,C_{2}$ in (5.2).

As the Shifted Ralston and Euler methods have the same orders of convergence
as the SRA1 and Euler-Maruyama methods, we estimate the ratios of their
$L^{2}(\mathbb{P})$ errors. Our results are presented in the following graph,
where we see the proposed methods achieve roughly a $3\times$ improvement in
their accuracy for sufficiently small step sizes.

Figure 5: $S_{N}$ estimated for (43) with 1,000,000 sample paths, where $N$ is the number of steps. To better illustrate differences in accuracy between methods, we plot the ratio $S_{N}^{(\text{Method 1})}/S_{N}^{(\text{Method 2)}}$. Table 3: Computer time to simulate $100,000$ sample paths of (43) in C++ with $100$ steps (seconds) Shifted Ralston (5.2) | SRA1 scheme [71] | Shifted Euler (42) | Euler-Maruyama
---|---|---|---
3.16 | 2.29 | 1.91 | 1.09

We highlight the $L^{2}(\mathbb{P})$ error ratio for $h=\frac{1}{1000}$ as it
can be explained theoretically. The shifted Ralston method uses the optimal
approximation (58), which has an error:

$\displaystyle\mathbb{E}\bigg{[}\bigg{(}\frac{1}{2}\int_{t_{k}}^{t_{k+1}}W_{t_{k}\hskip
0.7113pt,r}^{2}\,dr-\mathbb{E}\bigg{[}\frac{1}{2}\int_{t_{k}}^{t_{k+1}}W_{t_{k},r}^{2}\,dr\hskip
0.7113pt\Big{|}\hskip 0.7113ptW_{k}\hskip 0.7113pt,H_{k}\hskip
0.7113pt,n_{k}\bigg{]}\bigg{)}^{2}\,\bigg{]}^{\frac{1}{2}}$
$\displaystyle\hskip
14.22636pt=\mathbb{E}\bigg{[}\frac{11}{25200}h^{4}+\Big{(}\frac{1}{720}-\frac{1}{384\pi}\Big{)}h^{3}W_{k}^{2}+\frac{1}{700}h^{3}H_{k}^{2}-\frac{1}{320\sqrt{6\pi}}n_{k}h^{\frac{7}{2}}W_{k}\bigg{]}^{\frac{1}{2}}$
$\displaystyle\hskip
14.22636pt=\bigg{(}\frac{7}{3600}-\frac{1}{384\pi}\bigg{)}^{\frac{1}{2}}h^{2},$

where the second line follows from the conditional variance (59). On the other
hand, the Taylor expansion of the SRA1 scheme contains
$\frac{3}{16}(W_{k}+2H_{k})^{2}$, and has the error:

$\displaystyle\mathbb{E}\bigg{[}\bigg{(}\frac{1}{2}\int_{t_{k}}^{t_{k+1}}W_{t_{k},r}^{2}\,dr-\frac{3}{16}\big{(}W_{k}+2\hskip
0.7113ptH_{k}\big{)}^{2}\bigg{)}^{2}\,\bigg{]}^{\frac{1}{2}}$
$\displaystyle\hskip
14.22636pt=\mathbb{E}\bigg{[}\bigg{(}\frac{1}{48}hW_{k}^{2}+\frac{1}{4}hW_{k}H_{k}+\frac{3}{4}hH_{k}^{2}-L_{k}\bigg{)}^{2}\,\bigg{]}^{\frac{1}{2}}$
$\displaystyle\hskip
14.22636pt=\bigg{(}\mathbb{E}\bigg{[}\bigg{(}\frac{1}{48}hW_{k}^{2}+\cdots\bigg{)}^{2}\bigg{]}-2\hskip
0.7113pt\mathbb{E}\bigg{[}L_{k}\bigg{(}\frac{1}{48}hW_{k}^{2}+\frac{1}{4}hW_{k}H_{k}+\frac{3}{4}hH_{k}^{2}\bigg{)}\bigg{]}+\mathbb{E}\big{[}L_{k}^{2}\big{]}\bigg{)}^{\frac{1}{2}}$
$\displaystyle\hskip 14.22636pt=\bigg{(}\frac{1}{48}h^{4}-2\hskip
0.7113pt\mathbb{E}\bigg{[}\bigg{(}\frac{1}{30}h^{2}+\frac{3}{5}hH_{k}^{2}\bigg{)}\bigg{(}\frac{1}{48}hW_{k}^{2}+\frac{1}{4}hW_{k}H_{k}+\frac{3}{4}hH_{k}^{2}\bigg{)}\bigg{]}+\frac{1}{72}\hskip
0.7113pth^{4}\bigg{)}^{\frac{1}{2}}$ $\displaystyle\hskip
14.22636pt=\bigg{(}\frac{1}{120}\bigg{)}^{\frac{1}{2}}h^{2},$

where we used Theorems 3.9 and 3.10 from [29] to compute the above
expectations. Hence we can compute the ratio of these $L^{2}(\mathbb{P})$
errors as $\big{(}\frac{7}{30}-\frac{5}{16\pi}\big{)}^{\frac{1}{2}}=0.37$
(2.d.p). Unsurprisingly, this is close to the error ratio of 0.38 that was
seen in the experiment.

As both the Shifted Ralston (5.2) and SRA1 [71] methods require two
evaluations of the drift vector field per step, we expect them to have a
similar computational cost in settings where drift evaluations are expensive
(such as Langevin Monte Carlo [53]). However, if the noise is scalar, then the
scheme based on the non-linear splitting path (32) has the potential be almost
three times more accurate, for sufficiently small $h$.

More generally, when the Brownian motion is multidimensional, the non-linear
splitting approach is not able to accurately approximate the “cross” iterated
integrals

(44) $\displaystyle\underset{\begin{subarray}{c}\hskip 0.7113ptr_{3}\hskip
0.7113pt<\hskip 0.7113ptr_{2}\hskip 0.7113pt<\hskip 0.7113ptr_{1}\hskip
0.7113pt,\\\\[2.0pt] r_{i}\in[t_{k}\hskip
0.7113pt,t_{k+1}]\end{subarray}}{\int\int\int}dW_{r_{3}}^{\hskip
0.7113pti}\circ\hskip 0.7113ptdW_{r_{2}}^{\hskip
0.7113ptj}\,dr_{1}+\underset{\begin{subarray}{c}\hskip 0.7113ptr_{3}\hskip
0.7113pt<\hskip 0.7113ptr_{2}\hskip 0.7113pt<\hskip 0.7113ptr_{1}\hskip
0.7113pt,\\\\[2.0pt] r_{i}\in[t_{k}\hskip
0.7113pt,t_{k+1}]\end{subarray}}{\int\int\int}dW_{r_{3}}^{j}\circ\hskip
0.7113ptdW_{r_{2}}^{\hskip
0.7113pti}\,dr_{1}=\int_{t_{k}}^{t_{k+1}}W_{t_{k},r}^{\hskip
0.7113pti}W_{t_{k},r}^{\hskip 0.7113ptj}\hskip 0.7113ptdr,$

where $i\neq j$. On the other hand, we can see that the linear splitting path
(31) satisfies

$\mathbb{E}\bigg{[}\int_{t_{k}}^{t_{k+1}}\gamma_{t_{k},r}^{\hskip
0.7113pti}\gamma_{t_{k},r}^{\hskip 0.7113ptj}\hskip 0.7113ptdr\hskip
0.7113pt\Big{|}\hskip
0.7113ptW_{k}\bigg{]}=\mathbb{E}\bigg{[}\int_{t_{k}}^{t_{k+1}}W_{t_{k},r}^{\hskip
0.7113pti}W_{t_{k},r}^{\hskip 0.7113ptj}\hskip 0.7113ptdr\hskip
0.7113pt\Big{|}\hskip 0.7113ptW_{k}\bigg{]}=\frac{1}{3}hW_{k}^{\hskip
0.7113pti}W_{k}^{\hskip 0.7113ptj}+\frac{1}{6}h^{2}\delta_{ij}\hskip
0.7113pt,$

where the expectation is computable via standard properties of the Brownian
bridge. Numerical schemes for commutative SDEs matching the above conditional
expectation are known as “asymptotically efficient” (see [12, 15, 65] for
examples of such methods).

An attractive feature of both the shifted Euler and Ralston approaches is that
they are efficient in terms of drift evaluations per step, requiring one and
two respectively. This is particularly appealing in applications where drift
evaluations are expensive (such as in machine learning [45, 53]). We leave
further investigations and the analysis of “Shifted” Runge-Kutta methods for
additive noise SDEs as a topic of future work.

###  5.3 FitzHugh-Nagumo model

We consider a stochastic FitzHugh-Nagumo (FHN) model which has been used for
describing the spike activity of neurons [8, 51]. The stochastic FHN model
follows the two-dimensional additive noise SDE given by

(45) $d\begin{pmatrix}v_{t}\\\\[3.0pt]
u_{t}\end{pmatrix}=\begin{pmatrix}\frac{1}{\epsilon}\big{(}v_{t}-v_{t}^{3}-u_{t}\big{)}\\\\[3.0pt]
\gamma v_{t}-u_{t}+\beta\end{pmatrix}dt+\begin{pmatrix}\sigma_{1}&0\\\\[3.0pt]
0&\sigma_{2}\end{pmatrix}dW_{t}\hskip 0.7113pt.$

To discretise the stochastic FHN model, we apply the piecewise linear path
(28) and, similar to [8], apply a Strang splitting to approximate the
resulting drift ODE. Since $v^{\prime}=\frac{1}{\epsilon}(v-v^{3})$ admits a
closed-form solution, this leads to the splitting method:

$\displaystyle\begin{pmatrix}V_{k}^{(1)}\\\\[3.0pt] U_{k}^{(1)}\end{pmatrix}$
$\displaystyle:=\begin{pmatrix}V_{k}\\\\[3.0pt]
U_{k}\end{pmatrix}+\begin{pmatrix}\sigma_{1}&0\\\\[3.0pt]
0&\sigma_{2}\end{pmatrix}\begin{pmatrix}\frac{1}{2}W_{k}^{1}+H_{k}^{1}-\frac{1}{2}C_{k}^{1}\\\\[3.0pt]
\frac{1}{2}W_{k}^{2}+H_{k}^{2}-\frac{1}{2}C_{k}^{2}\end{pmatrix},$
$\displaystyle\begin{pmatrix}V_{k}^{(2)}\\\\[3.0pt] U_{k}^{(2)}\end{pmatrix}$
$\displaystyle:=\varphi_{\frac{1}{2}h}^{\text{Strang}}\begin{pmatrix}V_{k}^{(1)}\\\\[3.0pt]
U_{k}^{(1)}\end{pmatrix}+\begin{pmatrix}\sigma_{1}&0\\\\[3.0pt]
0&\sigma_{2}\end{pmatrix}\begin{pmatrix}C_{k}^{1}\\\\[3.0pt]
C_{k}^{2}\end{pmatrix},$ (46) $\displaystyle\begin{pmatrix}V_{k+1}\\\\[3.0pt]
U_{k+1}\end{pmatrix}$
$\displaystyle:=\varphi_{\frac{1}{2}h}^{\text{Strang}}\begin{pmatrix}V_{k}^{(2)}\\\\[3.0pt]
U_{k}^{(2)}\end{pmatrix}+\begin{pmatrix}\sigma_{1}&0\\\\[3.0pt]
0&\sigma_{2}\end{pmatrix}\begin{pmatrix}\frac{1}{2}W_{k}^{1}-H_{k}^{1}-\frac{1}{2}C_{k}^{1}\\\\[3.0pt]
\frac{1}{2}W_{k}^{2}-H_{k}^{2}-\frac{1}{2}C_{k}^{2}\end{pmatrix},$

where, for $u,v\in\mathbb{R}$, we define
$\varphi_{\frac{1}{2}h}^{\text{Strang}}\begin{pmatrix}v\\\ u\end{pmatrix}$ as

$\displaystyle\varphi_{\frac{1}{2}h}^{\text{Strang}}\begin{pmatrix}v\\\\[3.0pt]
u\end{pmatrix}$
$\displaystyle:=\begin{pmatrix}\widetilde{v}\Big{(}e^{-\frac{h}{2\epsilon}}+\widetilde{v}^{\hskip
0.7113pt2}\big{(}1-e^{-\frac{h}{2\epsilon}}\big{)}\Big{)}^{-\frac{1}{2}}\\\
\widetilde{u}+\frac{1}{4}\beta h\end{pmatrix},$

with $\widetilde{v}$ and $\widetilde{u}$ defined by

$\displaystyle\begin{pmatrix}\widetilde{v}\\\\[3.0pt]
\widetilde{u}\end{pmatrix}$
$\displaystyle:=\exp\Bigg{(}\frac{1}{2}h\begin{pmatrix}0&-\frac{1}{\epsilon}\\\\[3.0pt]
\gamma&-1\end{pmatrix}\Bigg{)}\begin{pmatrix}v\Big{(}e^{-\frac{h}{2\epsilon}}+v^{2}\big{(}1-e^{-\frac{h}{2\epsilon}}\big{)}\Big{)}^{-\frac{1}{2}}\\\
u+\frac{1}{4}\beta h\end{pmatrix},$

and the explicit formula for the above matrix exponential is given in [8,
Section 6.2]. The random variables $C_{k}^{1}$ and $C_{k}^{2}$ are given by
(29) and are generated independently. We note that, similar to the CIR model,
the stochastic FHN model is challenging to accurately simulate due to the
vector field not being globally Lipschitz continuous. That said, as the drift
does have polynomial growth and satisfies a one-sided Lipschitz condition,
there are numerical methods for (45) with strong convergence guarantees. We
will compare our scheme (5.3) against two such methods; the Strang splitting
scheme proposed in [8] and the Tamed Euler-Maruyama method introduced in [40].

Figure 6: $S_{N}$ estimated for (45) using 1,000 sample paths as a function of
step size $h=\frac{T}{N}$. The estimated strong errors for the Strang
splitting and Tamed Euler schemes were taken from [8].

In this numerical experiment, we used the following parameters in the FHN
model.

$\epsilon=1,\hskip 14.22636pt\gamma=1,\hskip 14.22636pt\beta=1,\hskip
14.22636pt\sigma_{1}=1,\hskip 14.22636pt\sigma_{2}=1,\hskip
14.22636pt(v_{0}\hskip 0.7113pt,u_{0})=(0,0),\hskip 14.22636ptT=5.$

We see in Figure 6 that the proposed high order splitting exhibits an 3/2
strong convergence rate and is significantly more accurate than other schemes
(for fixed $h$). For example, our splitting approach achieves better accuracy
in 320 steps than Strang splitting does in 10240 steps. As before, we present
simulation times for each method.

Table 4: Computer time to simulate $1,000$ sample paths of (45) in Python with $100$ steps (seconds) High order splitting (5.3) | Strang Splitting [8] | Tamed Euler-Maruyama [40]
---|---|---
8.15 | 2.66 | 1.71
Table 5: Estimated time to produce $1,000$ sample paths of (45) with an error of $S_{N}=10^{-3}$ (seconds) High order splitting (5.3) | Strang Splitting [8] | Tamed Euler-Maruyama [40]
---|---|---
10.4 | 110 | 166

From Tables 4 and 5, we conclude that the proposed high order splitting
method, which was derived using the path (28), gives the best performance for
the FHN model.

###  5.4 Underdamped Langevin dynamics

Traditionally used as a model for molecular dynamics [50, 59], the underdamped
Langevin diffusion (ULD) is given by

(47) $\displaystyle dx_{t}$ $\displaystyle=v_{t}\,dt,$ $\displaystyle dv_{t}$
$\displaystyle=-\gamma v_{t}\,dt-u\nabla f(x_{t})\,dt+\sqrt{2\gamma
u}\,dW_{t}\hskip 0.7113pt,$

where $x,v\in\mathbb{R}^{d}$ are the position and momentum of a particle,
$f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is a scalar potential, $\gamma>0$ is a
friction coefficient and $W$ is a $d$-dimensional Brownian motion. Under mild
conditions of $f$, the SDE (47) is known to admit a strong solution that is
ergodic with stationary distribution $\pi(x,v)\propto
e^{-f(x)}e^{-\frac{1}{2u}\|v\|^{2}}$ [69, Proposition 6.1]. As a consequence,
there has been recent interest in the application of (47) as an MCMC method
for high-dimensional sampling [9, 10, 13, 14, 19, 30, 39, 61, 72, 73, 76].

In this section, we briefly summarise the results of [30], where the path (34)
and a third order Runge-Kutta method are used to derive the following
numerical scheme:

$\displaystyle V_{n}^{(1)}$ $\displaystyle:=V_{n}+\sqrt{2\gamma
u}\,(H_{n}+6K_{n}),$ $\displaystyle X_{n}^{(1)}$
$\displaystyle:=X_{n}+\bigg{(}\frac{1-e^{-\frac{1}{2}\gamma
h}}{\gamma}\bigg{)}V_{n}^{(1)}-\bigg{(}\frac{e^{-\frac{1}{2}\gamma
h}+\frac{1}{2}\gamma h-1}{\gamma^{2}}\bigg{)}u\nabla f(X_{n})$
$\displaystyle\hskip 29.87538pt+\bigg{(}\frac{e^{-\frac{1}{2}\gamma
h}+\frac{1}{2}\gamma h-1}{\gamma^{2}h}\bigg{)}\sqrt{2\gamma
u}\,\big{(}W_{n}-12H_{n}\big{)},$ $\displaystyle X_{n+1}$
$\displaystyle:=X_{n}+\bigg{(}\frac{1-e^{-\gamma
h}}{\gamma}\bigg{)}V_{n}^{(1)}-\bigg{(}\frac{e^{-\gamma h}+\gamma
h-1}{\gamma^{2}}\bigg{)}\bigg{(}\frac{1}{3}u\nabla
f(X_{n})+\frac{2}{3}u\nabla\big{(}X_{n}^{(1)}\big{)}\bigg{)}$
$\displaystyle\hskip 29.87538pt+\bigg{(}\frac{e^{-\gamma h}+\gamma
h-1}{\gamma^{2}h}\bigg{)}\sqrt{2\gamma u}\,\big{(}W_{n}-12H_{n}\big{)},$
$\displaystyle V_{n}^{(2)}$ $\displaystyle:=e^{-\gamma
h}V_{n}^{(1)}-\frac{1}{6}e^{-\gamma h}u\nabla
f(X_{n})h-\frac{2}{3}e^{-\frac{1}{2}\gamma h}u\nabla
f\big{(}X_{n}^{(1)}\big{)}h-\frac{1}{6}u\nabla f(X_{n+1})h$
$\displaystyle\hskip 29.87538pt+\bigg{(}\frac{1-e^{-\gamma
h}}{\gamma}h\bigg{)}\sqrt{2\gamma u}\,\big{(}W_{n}-12H_{n}\big{)},$
$\displaystyle V_{n+1}$ $\displaystyle:=V_{n}^{(2)}-\sqrt{2\gamma
u}\,(H_{n}-6K_{n}),$

which we refer to as the SORT111Shifted ODE with Runge-Kutta Three method. It
is also worth noting that, since the final gradient evaluation $\nabla
f(X_{n+1})$ can be used in the next step to compute $(X_{n+2},V_{n+2})$, the
SORT method uses just two additional evaluations of the gradient $\nabla f$
per step. Moreover, as evaluating $\nabla f$ is usually much more
computationally expensive than generating $d$-dimensional Gaussian random
variables in practice, it follows that the SORT method is about twice as
expensive per step as the Euler-Maruyama method.

It was shown in [30] that under smoothness and convexity assumptions on $f$,
the Shifted ODE approximation based on (34) can achieve third order
convergence. However, this error analysis relies on properties of the ODE, and
thus obtaining such convergence guarantees for the SORT method itself remains
a topic of future research.

In the numerical experiment, we consider an application of ULD in data
science, namely the simulation of ULD as an MCMC algorithm for Bayesian
logistic regression. We use German credit data in [54], where each of the
$m=1000$ individuals has $d=49$ features $x_{i}\in\mathbb{R}^{d}$ and a label
$y_{i}\in\\{-1,1\\}$ indicating if they are creditworthy or not. The Bayesian
logistic regression model states that
$\mathbb{P}(Y_{i}=y_{i}|x_{i})=(1+e^{-y_{i}x_{i}^{\mathsf{T}}\theta})^{-1}$
where $\theta\in\mathbb{R}^{d}$ are parameters coming from the target density
$\pi(\theta)\propto\exp(-f(\theta))$ with

$f(\theta)=\delta\|\theta\|^{2}+\sum_{i=1}^{m}\log\big{(}1+\exp\big{(}-y_{i}x_{i}^{\mathsf{T}}\theta\big{)}\big{)}\hskip
0.7113pt.$

In the experiment, the regularisation parameter is set to $\delta=0.05$, the
ULD parameters are set to $u=1,\gamma=2$ and the initial condition
$\theta_{0}$ is sampled from a Gaussian prior,

$\theta_{0}\sim\mathcal{N}\big{(}0,10I_{d}\big{)}.$

In addition, we use a fixed time horizon of $T=1000$. The results are
presented below:

Figure 7: $S_{N}$ estimated for (47) with 1,000 sample paths as a function of
step size $h=\frac{T}{N}$.

From the above graph, we see that the SORT method [30] and UBU splitting [72]
are the best performing numerical schemes for ULD (in terms of strong
accuracy). We note that both the UBU [72] and Strang [9] splittings require a
single additional gradient evaluation per step and achieve essentially the
same accuracy in this example. Despite the SORT method having twice the
computational cost as UBU/Strang per step, we see that it is the best
performing method to achieve an accuracy of $S_{N}\leq 0.05$. Moreover, we
believe that the SORT method is the first approximation of ULD to exhibit
third order convergence whilst only requiring evaluations of the gradient
$\nabla f$.

##  6 Conclusion and future research

In this paper, we have presented a new simple methodology for designing and
analysing splitting methods for SDE simulation. The key idea is to replace the
system’s Brownian motion with a piecewise linear path. Moreover, for SDEs
satisfying a commutativity condition, we developed several high order
splitting methods which displayed state-of-the-art convergence in experiments.
As part of this investigation, we also detailed how recently developed
estimators for iterated integrals of Brownian motion can be directly
incorporated into such methods. Since these estimators were simply obtained as
the expectation of iterated integrals conditional on the generatable random
variables, they are optimal in an $L^{2}(\mathbb{P})$ sense. The technical
details underlying the integral estimators are presented in Appendix B.

Furthermore, the results in this paper may lead to several areas of future
research:

  * •

Development and analysis of methods inspired by splitting paths

For example, in the additive noise setting, the following Strong 1.5
Stochastic Runge-Kutta method is based on the splitting path (34), but with
$K_{s,t}=0$.

$\displaystyle\widetilde{Y}_{k}^{\text{SRK}}$
$\displaystyle:=Y_{k}^{\text{SRK}}+\sigma\hskip 0.7113ptH_{k}\hskip 0.7113pt,$
$\displaystyle\widetilde{Y}_{k+\frac{5}{6}}^{\text{SRK}}$
$\displaystyle:=\widetilde{Y}_{k}^{\text{SRK}}+\frac{5}{6}\Big{(}f\big{(}\widetilde{Y}_{k}^{\text{SRK}}\big{)}h+\sigma\hskip
0.7113ptW_{k}\Big{)}\hskip 0.7113pt,$ (48) $\displaystyle
Y_{k+1}^{\text{SRK}}$
$\displaystyle:=Y_{k}^{\text{SRK}}+\frac{2}{5}f\big{(}\widetilde{Y}_{k}^{\text{SRK}}\big{)}h+\frac{3}{5}f\big{(}\widetilde{Y}_{k+\frac{5}{6}}^{\text{SRK}}\big{)}h+\sigma\hskip
0.7113ptW_{k}\hskip 0.7113pt.$

As (• ‣ 6) does not follow a high order splitting path or use Ralston’s
method, it was not included in Section 5.2 and thus we leave its analysis as
future work. Similarly, conducting error analyses and further numerical
investigations for the Shifted Euler and Runge-Kutta methods from Section 5 is
a future topic.

  * •

Development of high order splitting methods for general SDEs

For example, the below method is a combination of the Strang splitting (26)
and the log-ODE method from rough path theory [63, Appendices A and B].

$\displaystyle Y_{k+1}:=\exp\bigg{(}\frac{1}{2}\hskip 0.7113ptf(\hskip
0.7113pt\cdot\hskip 0.7113pt)\hskip 0.7113pth\bigg{)}\exp\bigg{(}g(\hskip
0.7113pt\cdot\hskip 0.7113pt)W_{k}+\sum_{i\hskip 0.7113pt<\hskip
0.7113ptj}\big{[}\hskip 0.7113ptg_{i}\hskip 0.7113pt,g_{j}\big{]}(\hskip
0.7113pt\cdot\hskip 0.7113pt)A_{k}^{ij}\bigg{)}\exp\bigg{(}\frac{1}{2}\hskip
0.7113ptf(\hskip 0.7113pt\cdot\hskip 0.7113pt)\hskip
0.7113pth\bigg{)}Y_{k}\hskip 0.7113pt,$

where $[\hskip 0.7113ptg_{i}\hskip
0.7113pt,g_{j}](\cdot)=g_{j}^{\prime}(\cdot)g_{i}(\cdot)-g_{i}^{\prime}(\cdot)g_{j}(\cdot)$
is the standard vector field Lie bracket and $A_{k}=\\{A_{k}^{ij}\\}_{1\leq
i,j\leq d}$ is the Lévy area of the Brownian motion over $[t_{k}\hskip
0.7113pt,t_{k+1}]$. If $A_{k}$ is replaced by a random matrix
$\widetilde{A}_{k}\hskip 0.7113pt$, with the same mean and covariance, we
expect the resulting splitting method to achieve $O(h^{2})$ weak convergence.
Similarly, we expect that the Ninomiya-Ninomiya [66] and Ninomiya-Victoir [67]
weak second order schemes can be reinterpreted as path-based splittings.

  * •

Incorporating $\boldsymbol{(W_{k}\hskip 0.7113pt,H_{k}\hskip
0.7113pt,n_{k})}$-methods into Multilevel Monte Carlo

Multilevel Monte Carlo (MLMC), introduced by Giles in [35], is a popular
control variate strategy for variance reduction in SDE simulation. Whilst the
MLMC literature is quite extensive, we refer the reader to [1, 21, 34, 36] for
details on some of the improvements to MLMC that have been proposed. The
variance reduction obtained by MLMC is a consequence of the strong convergence
properties of the SDE solver as sample paths are computed for each level using
two step sizes ($h$ and $\frac{1}{2}h$), but with the same Brownian paths.
However, the majority of MLMC methods for SDE simulation use only the
increments of the Brownian motion. Thus, we conjecture that further variance
reduction can be achieved using high order $(W_{k}\hskip 0.7113pt,H_{k}\hskip
0.7113pt,n_{k})$-based SDE solvers. More generally, the application of high
order numerical methods to MLMC has seen relatively little attention (see [1,
21] for examples of such research).

  * •

Incorporating adaptive step sizes into $\boldsymbol{(W_{k}\hskip
0.7113pt,H_{k}\hskip 0.7113pt,n_{k})}$-based methods

Since it is possible to generate both $(W_{s,u}\hskip 0.7113pt,H_{s,u}\hskip
0.7113pt,n_{s,u})$ and $(W_{u,t}\hskip 0.7113pt,H_{u,t}\hskip
0.7113pt,n_{u,t})$ conditional on $(W_{s,t}\hskip 0.7113pt,H_{s,t}\hskip
0.7113pt,n_{s,t})$, where $u=s+\frac{1}{2}h$ is the midpoint of $[s,t]$, the
proposed splitting methods can be applied using an adaptive step size. Such a
methodology was detailed and initial investigated in [26, Chapter 6].

  * •

Application to high-dimensional SDEs in physics and data science

High-dimensional SDEs have seen a variety of real-world applications, ranging
from molecular dynamics [50, 59] to machine learning [44, 45, 52, 75, 79, 82].
Therefore, it would be interesting to investigate whether the splitting
methods developed in this paper could improve algorithms used in these
applications.

## References

  * [1] A. Al Gerbi, B. Jourdain, and E. Clément, Ninomiya–Victoir scheme: Strong convergence, antithetic version and application to multilevel estimators, Monte Carlo Methods and Applications, 22 (2016), pp. 197–228.
  * [2] A. Alfonsi, On the discretization schemes for the CIR (and Bessel squared) processes, Monte Carlo Methods and Applications, 11 (2005), pp. 355–384.
  * [3] A. Alfonsi, High order discretization schemes for the CIR process: Application to Affine Term Structure and Heston models, Mathematics of Computation, 79 (2010), pp. 209–237.
  * [4] A. Alfonsi, Strong order one convergence of a drift implicit Euler scheme: Application to the CIR process, Statistics & Probability Letters, 83 (2013), pp. 602–607.
  * [5] C. Bayer, The geometry of iterated stratonovich integrals, preprint, (2006), https://www.wias-berlin.de/people/bayerc/files/strat_geom.pdf.
  * [6] S. Biswas, C. Kumar, Neelima, G. d. Reis, and C. Reisinger, An explicit Milstein-type scheme for interacting particle systems and McKean–Vlasov SDEs with common noise and non-differentiable drift coefficients, https://arxiv.org/2208.10052, (2022).
  * [7] D. Brigo and F. Mercurio, Interest Rate Models – Theory and Practice, second edition, Springer, 2006.
  * [8] E. Buckwar, A. Samson, M. Tamborrino, and I. Tubikanec, A splitting method for SDEs with locally Lipschitz drift: Illustration on the FitzHugh-Nagumo model, Applied Numerical Mathematics, 179 (2022), pp. 191–220.
  * [9] E. Buckwar, M. Tamborrino, and I. Tubikanec, Spectral density-based and measure-preserving ABC for partially observed diffusion processes. An illustration on Hamiltonian SDEs, Statistics and Computing, 30 (2020), pp. 627–648.
  * [10] G. Bussi and M. Parrinello, Accurate sampling using Langevin dynamics, Physical Review E, 75 (2007).
  * [11] T. Cass and P. Friz, Densities for rough differential equations under hörmander’s condition, Annals of mathematics, (2010), pp. 2115–2141.
  * [12] F. Castell and J. Gaines, The ordinary differential equation approach to asymptotically efficient schemes for solution of stochastic differential equations, Annales de l’Institut Henri Poincaré, Probabilités et Statistiques, 32 (1996), pp. 231–250.
  * [13] M. Chak, N. Kantas, T. Lelièvre, and G. A. Pavliotis, Optimal friction matrix for underdamped Langevin sampling, https://arxiv.org/abs/2112.06844, (2021).
  * [14] X. Cheng, N. S. Chatterji, P. L. Bartlett, and M. I. Jordan, Underdamped Langevin MCMC: A non-asymptotic analysis, Proceedings of the 31st Conference On Learning Theory, Volume 75 of Proceedings of Machine Learning Research, (2018).
  * [15] J. M. C. Clark, An efficient approximation scheme for a class of stochastic differential equations, in Advances in Filtering and Optimal Stochastic Control, Springer, 1982.
  * [16] J. M. C. Clark and R. J. Cameron, The maximum rate of convergence of discrete approximations for stochastic differential equations, in Stochastic Differential Systems Filtering and Control, ed. by Grigelionis (Springer, Berlin), 1980.
  * [17] J. C. Cox, J. E. Ingersoll, and S. A. Ross, A Theory of the Term Structure of Interest Rates, Econometrica, 53 (1985), pp. 385–407.
  * [18] A. Cozma and C. Reisinger, Strong order 1/2 convergence of full truncation Euler approximations to the Cox–Ingersoll–Ross process, IMA Journal of Numerical Analysis, 40 (2020), pp. 358–376.
  * [19] A. S. Dalalyan and L. Riou-Durand, On sampling from a log-concave density using kinetic Langevin diffusions, Bernoulli, 26 (2020), pp. 1956–1988.
  * [20] A. Davie, KMT theory applied to approximations of SDE, in Stochastic Analysis and Applications, vol. 100 of Springer Proceedings in Mathematics and Statistics, Springer, 2014, pp. 185–201.
  * [21] K. Debrabant and A. Rößler, On the Acceleration of the Multi-Level Monte Carlo Method, Journal of Applied Probability, 52 (2018), pp. 307–322.
  * [22] S. Dereich, A. Neuenkirch, and L. Szpruch, An Euler-type method for the strong approximation of the Cox–Ingersoll–Ross process, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 468 (2011), pp. 1105–1115.
  * [23] A. S. Dickinson, Optimal Approximation of the Second Iterated Integral of Brownian Motion, Stochastic Analysis and Applications, 25 (2007), pp. 1109–1128.
  * [24] R. C. Elandt, The Folded Normal Distribution: Two Methods of Estimating Parameters from Moments, Technometrics, 3 (1961), pp. 551–562.
  * [25] W. Fang and M. B. Giles, Adaptive Euler-Maruyama method for SDEs with nonglobally Lipschitz drift, Annals of Applied Probability, 30 (2020), pp. 526–560.
  * [26] J. Foster, Numerical approximations for stochastic differential equations, PhD thesis, University of Oxford, 2020, https://ora.ox.ac.uk/objects/uuid:775fc3f5-501c-425f-8b43-fc5a7b2e4310.
  * [27] J. Foster and K. Habermann, Brownian bridge expansions for Lévy area approximations and particular values of the Riemann zeta function, Combinatorics, Probability and
Computing, (2022).

  * [28] J. Foster, T. Lyons, and V. Margarint, An asymptotic radius of convergence for the Loewner equation and simulation of SLE traces via splitting, Journal of Statistical Physics, 189 (2022).
  * [29] J. Foster, T. Lyons, and H. Oberhauser, An Optimal Polynomial Approximation of Brownian Motion, SIAM Journal on Numerical Analysis, 58 (2020), pp. 1393–1421.
  * [30] J. Foster, T. Lyons, and H. Oberhauser, The shifted ODE method for underdamped Langevin MCMC, https://arxiv.org/abs/2101.03446, (2021).
  * [31] P. K. Friz and M. Hairer, A Course on Rough Paths: With an Introduction to Regularity Structures, Springer, 2020.
  * [32] P. K. Friz and N. B. Victoir, Multidimensional Stochastic Processes as Rough Paths: Theory and Applications, vol. 120, Cambridge University Press, 2010.
  * [33] J. G. Gaines and T. Lyons, Random Generation of Stochastic Area Integrals, SIAM Journal on Applied Mathematics, 54 (1994), pp. 1132–1146.
  * [34] M. B. Giles, Improved multilevel Monte Carlo convergence using the Milstein scheme, in Monte Carlo and Quasi-Monte Carlo Methods, ed. by Keller, Heinrich and Niederreiter (Springer, Berlin), 2008.
  * [35] M. B. Giles, Multilevel Monte Carlo path simulation, Operations Research, 56 (2008), pp. 607–617.
  * [36] M. B. Giles and L. Szpruch, Antithetic multilevel Monte Carlo estimation for multi-dimensional SDEs without Lévy area simulation, Annals of Applied Probability, 24 (2014), pp. 1585–1620.
  * [37] M. Hefter and A. Jentzen, On arbitrarily slow convergence rates for strong numerical approximations of Cox-Ingersoll-Ross processes and squared Bessel processes, Finance and Stochastics, 23 (2019), pp. 139–172.
  * [38] S. L. Heston, A Closed-Form Solution for Options with Stochastic Volatility with Applications to Bond and Currency Options, The Review of Financial Studies, 6 (1993), pp. 327–343.
  * [39] Z. Hu, F. Huang, and H. Huang, Optimal Underdamped Langevin MCMC Method, Advances in Neural Information Processing Systems, (2021).
  * [40] M. Hutzenthaler, A. Jentzen, and P. E. Kloeden, Strong convergence of an explicit numerical method for SDEs with nonglobally Lipschitz continuous coefficients, Annals of Applied Probability, 22 (2012), pp. 1611–1641.
  * [41] Y. Iguchi and T. Yamada, Operator splitting around Euler-Maruyama scheme and high order discretization of heat kernels, ESAIM: Mathematical Modelling and Numerical Analysis, 55 (2021), pp. 323–367.
  * [42] C. Kelly, G. Lord, and H. Maulana, The role of adaptivity in a numerical method for the Cox–Ingersoll–Ross model, Journal of Computational and Applied Mathematics, 410 (2022).
  * [43] P. Kidger, On Neural Differential Equations, PhD thesis, University of Oxford, 2021, https://ora.ox.ac.uk/objects/uuid:af32d844-df84-4fdc-824d-44bebc3d7aa9.
  * [44] P. Kidger, J. Foster, X. Li, and T. Lyons, Efficient and Accurate Gradients for Neural SDEs, Advances in Neural Information Processing Systems, (2021).
  * [45] P. Kidger, J. Foster, X. Li, H. Oberhauser, and T. Lyons, Neural SDEs as Infinite-Dimensional GANs, Proceedings of the 38th International Conference on Machine Learning, (2021).
  * [46] P. Kidger, J. Morrill, J. Foster, and T. Lyons, Neural Controlled Differential Equations for Irregular Time Series, Advances in Neural Information Processing Systems, (2020).
  * [47] P. E. Kloeden and E. Platen, Numerical Solution of Stochastic Differential Equations, Springer, Berlin, 1992.
  * [48] R. Kruse and Y. Wu, A randomized and fully discrete Galerkin finite element method for semilinear stochastic evolution equations, Mathematics of Computation, 88 (2019), pp. 2793–2825.
  * [49] R. Kruse and Y. Wu, A randomized Milstein method for stochastic differential equations with non-differentiable drift coefficients, Discrete and Continuous Dynamical Systems. Series B. A Journal Bridging Mathematics and Sciences, 24 (2019), pp. 3475–3502.
  * [50] B. Leimkuhler and C. Matthews, Molecular Dynamics: With Deterministic and Stochastic Numerical Methods, Interdisciplinary Applied Mathematics, Springer, 2015.
  * [51] J. R. León and A. Samson, Hypoelliptic stochastic FitzHugh-Nagumo neuronal model: mixing, up-crossing and estimation of the spike rate, Annals of Applied Probability, 28 (2018), pp. 2243–2274.
  * [52] R. Li, H. Zha, and M. Tao, Sqrt(d) dimension dependence of Langevin Monte Carlo, Proceedings of the 10th International Conference on Learning Representations, (2022).
  * [53] X. Li, D. Wu, L. Mackey, and M. A. Erdogdu, Stochastic Runge-Kutta Accelerates Langevin Monte Carlo and Beyond, Advances in Neural Information Processing Systems, (2019).
  * [54] M. Lichman, UCI machine learning repository, https://archive.ics.uci.edu/ml, 2013.
  * [55] T. Lyons, M. Caruana, and T. Lévy, Differential Equations Driven by Rough Paths, vol. 1908 of Lecture Notes in Mathematics., Springer, 2007.
  * [56] T. Lyons and N. Victoir, Cubature on Wiener space, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 460 (2004), pp. 169–198.
  * [57] X. Mao, Stochastic Differential Equations and Applications, second edition, Elsevier, 2008.
  * [58] G. N. Milstein and J. Schoenmakers, Uniform approximation of the Cox-Ingersoll-Ross process, Advances in Applied Probability, 47 (2015), pp. 1132–1156.
  * [59] G. N. Milstein and M. V. Tretyakov, Stochastic Numerics for Mathematical Physics, Springer, 2004.
  * [60] T. Misawa, Numerical integration of stochastic differential equations by composition methods, in Dynamical systems and differential geometry (Japanese), no. 1180, 2000, pp. 166–190.
  * [61] P. Monmarché, High-dimensional MCMC with a standard splitting scheme for the underdamped Langevin diffusion, Electronic Journal of Statistics, 15 (2021), pp. 4117–4166.
  * [62] J. Morrill, P. Kidger, L. Yang, and T. Lyons, On the Choice of Interpolation Scheme for Neural CDEs, Transactions on Machine Learning Research, (2022).
  * [63] J. Morrill, C. Salvi, P. Kidger, J. Foster, and T. Lyons, Neural Rough Differential Equations for Long Time Series, Proceedings of the 38th International Conference on Machine Learning, (2021).
  * [64] J. Mrongowius and A. Rößler, On the approximation and simulation of iterated stochastic integrals and the corresponding Lévy areas in terms of a multidimensional Brownian motion, Stochastic Analysis and Applications, 40 (2022), pp. 397–425.
  * [65] N. J. Newton, Asymptotically Efficient Runge-Kutta Methods for a Class of Itô and Stratonovich Equations, SIAM Journal on Applied Mathematics, 51 (1991), pp. 303–604.
  * [66] M. Ninomiya and S. Ninomiya, A new higher-order weak approximation scheme for stochastic differential equations and the Runge–Kutta method, Finance and Stochastics, 13 (2009), pp. 415–443.
  * [67] S. Ninomiya and N. Victoir, Weak Approximation of Stochastic Differential Equations and Application to Derivative Pricing, Applied Mathematical Finance, 15 (2008), pp. 107–121.
  * [68] S. Pathiraja, L2 convergence of smooth approximations of Stochastic Differential Equations with unbounded coefficients, https://arxiv.org/abs/2011.13009, (2020).
  * [69] G. A. Pavliotis, Stochastic Processes and Applications, Springer, New York, 2014.
  * [70] A. Ralston, Runge-Kutta methods with minimum error bounds, Mathematics of Computation, (1962), pp. 431–437.
  * [71] A. Rößler, Runge–Kutta methods for the strong approximation of solutions of stochastic differential equations, SIAM Journal on Numerical Analysis, 8 (2010), pp. 922–952.
  * [72] J. M. Sanz-Serna and K. C. Zygalakis, Wasserstein distance estimates for the distributions of numerical approximations to ergodic stochastic differential equations, Journal of Machine Learning Research, 22 (2021).
  * [73] R. Shen and Y. T. Lee, The Randomized Midpoint Method for Log-Concave Sampling, Advances in Neural Information Processing Systems, (2019).
  * [74] A. Shmatkov, Rate of Convergence of Wong-Zakai Approximations for SDEs and SPDEs, PhD thesis, University of Edinburgh, 2005.
  * [75] Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, and B. Poole, Score-Based Generative Modeling through Stochastic Differential Equations, Proceedings of the International Conference on Learning Representations, (2021).
  * [76] Z. Song and Z. Tan, Hamiltonian-Assisted Metropolis Sampling, Journal of the American Statistical Association, Theory and Methods, (2021).
  * [77] G. Strang, On the construction and comparison of difference schemes, SIAM Journal on Numerical Analysis, 5 (1968), pp. 506–517.
  * [78] I. Tubikanec, M. Tamborrino, P. Lansky, and E. Buckwar, Qualitative properties of different numerical methods for the inhomogeneous geometric Brownian motion, Journal of Computational and Applied Mathematics, 406 (2022).
  * [79] M. Welling and Y. W. Teh, Bayesian Learning via Stochastic Gradient Langevin Dynamics, Proceedings of the 28th International Conference on Machine Learning, (2011).
  * [80] M. Wiktorsson, Joint characteristic function and simultaneous simulation of iterated Itô integrals for multiple independent Brownian motions, Annals of Applied Probability, 11 (2001), pp. 470–487.
  * [81] E. Wong and M. Zakai, On the Convergence of Ordinary Integrals to Stochastic Integrals, Annals of Mathematical Statistics, 36 (1965), pp. 1560–1564.
  * [82] Q. Zhang and Y. Chen, Path Integral Sampler: A Stochastic Control Approach For Sampling, Proceedings of the International Conference on Learning Representations, (2022).

##  Appendix A Cancellation of integrals under commutativity condition

In this section, we present the calculations required to simplify the Taylor
expansions of the SDE (1) and CDE (3) when the commutativity condition (2) is
satisfied. To make these calculations easier, we identify iterated integrals
with non-commutative polynomials and present the shuffle product, commonly
used in rough path theory [55].

######  Definition A.1.

Let $\mathcal{A}_{d}$ denote the set of letters $\\{0,1,\cdots,d\\}$. We can
identify linear combinations of iterated integrals with elements in
$\mathbb{R}\langle\mathcal{A}_{d}\rangle$ by $I_{e}^{W}=I_{e}^{\gamma}=1$ and

$\displaystyle i_{1}\cdots\hskip 0.7113pti_{m}$ $\displaystyle\leftrightarrow
I_{i_{1}\cdots\hskip
0.7113pti_{m}}^{W}:=\int_{s}^{t}\int_{s}^{r_{1}}\dots\int_{s}^{r_{m-1}}\circ\,dW^{i_{1}}_{r_{n}}\circ
dW^{i_{2}}_{r_{m-1}}\cdots\circ dW^{i_{m-1}}_{r_{2}}\circ
dW^{i_{m}}_{r_{1}}\hskip 0.7113pt,$ $\displaystyle i_{1}\cdots\hskip
0.7113pti_{m}$ $\displaystyle\leftrightarrow I_{i_{1}\cdots\hskip
0.7113pti_{m}}^{\gamma}:=\int_{s}^{t}\int_{s}^{r_{1}}\dots\int_{s}^{r_{m-1}}d\gamma^{i_{1}}_{r_{m}}\hskip
0.7113ptd\gamma^{i_{2}}_{r_{m-1}}\cdots\hskip
0.7113ptd\gamma^{i_{m-1}}_{r_{2}}\hskip 0.7113ptd\gamma^{i_{m}}_{r_{1}}\hskip
0.7113pt,$ $\displaystyle\lambda_{1}u+\lambda_{2}v$
$\displaystyle\leftrightarrow
I_{\lambda_{1}u+\lambda_{2}v}:=\lambda_{1}I_{u}+\lambda_{2}I_{v}\hskip
0.7113pt,$

for $m\geq 0$, $i_{1}\hskip 0.7113pt,i_{2}\hskip
0.7113pt,\cdots,i_{m}\in\mathcal{A}_{d}\hskip 0.7113pt$,
$u,v\in\mathcal{A}_{d}^{\ast}$ and $\lambda_{1}\hskip
0.7113pt,\lambda_{2}\in\mathbb{R}$.

######  Definition A.2.

Suppose that $\mathcal{A}_{d}$ is a set containing $d$ letters and let
$\mathbb{R}\langle\mathcal{A}_{d}\rangle$ be the corresponding space of non-
commutative polynomials in $\mathcal{A}_{d}$ with real coefficients. Then the
shuffle product $\hskip
0.7113pt\sqcup\mathchoice{\mkern-3.0mu}{\mkern-3.0mu}{\mkern-1.0mu}{\mkern-1.0mu}\sqcup:\mathbb{R}\langle\mathcal{A}_{d}\rangle\times\mathbb{R}\langle\mathcal{A}_{d}\rangle\rightarrow\mathbb{R}\langle\mathcal{A}_{d}\rangle$
is the unique bilinear map such that

$\displaystyle
ua\sqcup\mathchoice{\mkern-7.0mu}{\mkern-7.0mu}{\mkern-3.5mu}{\mkern-3.5mu}\sqcup\hskip
2.13394pt\,vb$
$\displaystyle=(u\sqcup\mathchoice{\mkern-7.0mu}{\mkern-7.0mu}{\mkern-3.5mu}{\mkern-3.5mu}\sqcup\hskip
2.13394ptvb)\hskip
0.7113pta+(ua\sqcup\mathchoice{\mkern-7.0mu}{\mkern-7.0mu}{\mkern-3.5mu}{\mkern-3.5mu}\sqcup\hskip
2.13394ptv)\hskip 0.7113ptb,$ $\displaystyle
u\sqcup\mathchoice{\mkern-7.0mu}{\mkern-7.0mu}{\mkern-3.5mu}{\mkern-3.5mu}\sqcup\hskip
2.13394pte$
$\displaystyle=e\sqcup\mathchoice{\mkern-7.0mu}{\mkern-7.0mu}{\mkern-3.5mu}{\mkern-3.5mu}\sqcup\hskip
2.13394ptu=u,$

where $e$ denotes the empty letter.

With this notation, we can link the shuffle project to the integration by
parts formula. As a result, the shuffle project will allow us to expand
products of iterated integrals.

######  Theorem A.3 (Integration by parts formula for integrals).

For all $\hskip 0.7113ptu,v\in\mathbb{R}\langle\mathcal{A}_{d}\rangle$, we
have

(49) $I_{u}\cdot I_{v}=I_{u\hskip 0.7113pt\sqcup\sqcup\hskip 0.7113ptv}$

######  Proof A.4.

It is clear that the identity (49) holds when $u=e$ or $v=e$ since $I_{e}=1$.
Suppose that (49) holds for all words $u,v\in\mathcal{A}_{d}^{\ast}$ with a
combined length less than $m$. Then for words
$u,v\in\mathcal{A}_{d}^{\ast}\hskip 0.7113pt$ and letters
$a,b\in\mathcal{A}_{d}\hskip 0.7113pt$ such that $|ua|+|vb|=m$, we have

$\displaystyle I_{ua}^{W}\cdot I_{vb}^{W}$
$\displaystyle=\int_{s}^{t}I_{u}^{W}(r)\circ
dW_{r}^{a}\int_{s}^{t}I_{v}^{W}(r)\circ dW_{r}^{b}$
$\displaystyle=\int_{s}^{t}\bigg{(}\int_{s}^{r_{1}}I_{u}^{W}(r_{2})\circ
dW_{r_{2}}^{a}\bigg{)}\circ d\bigg{(}\int_{s}^{r_{1}}I_{v}^{W}(r_{2})\circ
dW_{r_{2}}^{b}\bigg{)}$ $\displaystyle\hskip
14.22636pt+\int_{s}^{t}\bigg{(}\int_{s}^{r_{1}}I_{v}^{W}(r_{2})\circ
dW_{r_{2}}^{b}\bigg{)}\circ d\bigg{(}\int_{s}^{r_{1}}I_{u}^{W}(r_{2})\circ
dW_{r_{2}}^{a}\bigg{)}$ $\displaystyle=\int_{s}^{t}I_{ua}^{W}(r_{1})\hskip
0.7113ptI_{v}^{W}(r_{1})\circ
dW_{r_{1}}^{b}+\int_{s}^{t}I_{vb}^{W}(r_{1})\hskip
0.7113ptI_{u}^{W}(r_{1})\circ dW_{r_{1}}^{a}$ $\displaystyle=I_{(ua\hskip
0.7113pt\sqcup\sqcup\hskip 0.7113ptv)\hskip 0.7113ptb+(u\hskip
0.7113pt\sqcup\sqcup\hskip 0.7113ptvb)\hskip 0.7113pta}^{W}\hskip 0.7113pt,$

where the second line uses integration by parts (which holds for Stratonovich
integrals) and the last line uses the induction hypothesis. The result now
follows by linearity. The same argument gives (49) for iterated integrals with
respect to the path $\gamma\hskip 0.7113pt$.

Using Theorem A.3, it will be straightforward to rewrite products of integrals
as linear combinations of (high order) integrals. In addition, it shall enable
us to establish decompositions of iterated integrals into symmetric and
antisymmetric components.

######  Theorem A.5 (Symmetric and antisymmetric components of iterated
integrals).

Let the Lie bracket $\hskip 0.7113pt[\,\cdot\hskip 0.7113pt,\hskip
0.7113pt\cdot\hskip
0.7113pt]:\mathbb{R}\langle\mathcal{A}_{d}\rangle\times\mathbb{R}\langle\mathcal{A}_{d}\rangle\rightarrow\mathbb{R}\langle\mathcal{A}_{d}\rangle$
be the unique bilinear map with

(50) $[u,v]=uv-vu\hskip 0.7113pt,$

for words $u,v\in\mathcal{A}_{d}^{\ast}$. Then, adopting the notation of
Definition A.1 and Theorem A.3, we have

(51) $\displaystyle I_{ij}$ $\displaystyle=\frac{1}{2}I_{i}\cdot
I_{j}+\frac{1}{2}I_{[i,j]}\hskip 0.7113pt,$ (52) $\displaystyle I_{ijk}$
$\displaystyle=\frac{1}{6}I_{i}\cdot I_{j}\cdot I_{k}+\frac{1}{4}I_{i}\cdot
I_{[j,k]}+\frac{1}{4}I_{[i,j]}\cdot
I_{k}+\frac{1}{6}I_{[[i,j],k]}+\frac{1}{6}I_{[i,[j,k]]}\hskip 0.7113pt,$ (53)
$\displaystyle I_{ijkl}$ $\displaystyle=\frac{1}{24}I_{i}\cdot I_{j}\cdot
I_{k}\cdot I_{l}+\frac{1}{12}I_{i}\cdot I_{[j,[k,l]]}+\frac{1}{12}I_{i}\cdot
I_{[[j,k],l]}+\frac{1}{12}I_{[i,[j,k]]}\cdot I_{l}$ $\displaystyle\hskip
14.22636pt+\frac{1}{12}I_{[[i,j],k]}\cdot I_{l}+\frac{1}{12}I_{i}\cdot
I_{j}\cdot I_{[k,l]}+\frac{1}{12}I_{i}\cdot I_{[j,k]}\cdot
I_{l}+\frac{1}{12}I_{[i,j]}\cdot I_{k}\cdot I_{l}$ $\displaystyle\hskip
14.22636pt+\frac{1}{8}I_{[i,j]}\cdot
I_{[k,l]}+\frac{1}{12}I_{[i,[j,[k,l]]]}+\frac{1}{12}I_{[[i,[j,k]],l]}+\frac{1}{12}I_{[[[i,j],k],l]}$
$\displaystyle\hskip
14.22636pt+\frac{1}{12}\big{(}I_{kjli}-I_{kjil}+I_{lijk}-I_{iljk}\big{)}+\frac{1}{12}\big{(}I_{jilk}-I_{kilj}+I_{jlik}-I_{klij}\big{)}\hskip
0.7113pt,$

for $i,j,k,l\in\mathcal{A}_{d}\hskip 0.7113pt$.

######  Proof A.6.

The results follow by expanding the Lie brackets $[\,\cdot\hskip
0.7113pt,\hskip 0.7113pt\cdot\hskip 0.7113pt]$ in the right hand sides using
(50) and applying the integration by parts formula via Theorem 49.

Hence, in order to simplify the Taylor expansions of (1) and (3), we will need
to find symmetries in the vector field derivatives that will cause the
antisymmetric parts in the iterated integrals to cancel out. To this end, we
give the following lemma, which will set up the notation used in the main
result of the section (Theorem A.9).

######  Lemma A.7 (Codomains of vector field derivatives).

Given a sufficiently smooth vector field $\hskip
0.7113ptf:\mathbb{R}^{e}\rightarrow\mathbb{R}^{e}$, its Fréchet derivatives
will map between the following spaces:

$\displaystyle f^{\hskip 0.7113pt\prime}$
$\displaystyle:\mathbb{R}^{e}\rightarrow L(\mathbb{R}^{e},\mathbb{R}^{e}),$
$\displaystyle f^{\hskip 0.7113pt\prime\prime}$
$\displaystyle:\mathbb{R}^{e}\rightarrow
L(\mathbb{R}^{e},L(\mathbb{R}^{e},\mathbb{R}^{e})),$ $\displaystyle f^{\hskip
0.7113pt\prime\prime\prime}$ $\displaystyle:\mathbb{R}^{e}\rightarrow
L(\mathbb{R}^{e},L(\mathbb{R}^{e},L(\mathbb{R}^{e},\mathbb{R}^{e}))),$

where $L(U,V)$ denotes the space of linear maps between the vector spaces $U$
and $V$. Equivalently, we can view $f^{(k)}(y)$ as a $k$-linear map on
$\mathbb{R}^{e}$ for each $y\in\mathbb{R}^{e}$. That is,

$\displaystyle f^{\hskip 0.7113pt\prime}$
$\displaystyle:\mathbb{R}^{e}\rightarrow L(\mathbb{R}^{e},\mathbb{R}^{e}),$
$\displaystyle f^{\hskip 0.7113pt\prime\prime}$
$\displaystyle:\mathbb{R}^{e}\rightarrow L\big{(}(\mathbb{R}^{e})^{\otimes
2},\mathbb{R}^{e}\big{)},$ $\displaystyle f^{\hskip
0.7113pt\prime\prime\prime}$ $\displaystyle:\mathbb{R}^{e}\rightarrow
L\big{(}(\mathbb{R}^{e})^{\otimes 3},\mathbb{R}^{e}\big{)}.$

######  Proof A.8.

The result follows immediately from the definition of Fréchet derivative.

We now turn to the main result of the section, which shows that certain terms
in the Taylor expansions of (1) and (3) simplify under the commutativity
condition. We note that the CDE (3) can be written in the equivalent form

$\displaystyle\hskip 7.11317ptdy_{r}^{\gamma}=f(y_{r}^{\gamma})\hskip
0.7113ptd\gamma^{\tau}(r)+\sum_{i=1}^{d}g_{i}(y_{r}^{\gamma})\hskip
0.7113ptd(\gamma^{\hskip 0.7113pt\omega}(r))^{i}.$

######  Theorem A.9.

Suppose that the following commutativity condition holds

(54) $\hskip 7.11317ptg_{i}^{\hskip 0.7113pt\prime}(y)g_{j}(y)=g_{j}^{\hskip
0.7113pt\prime}(y)g_{i}(y),\hskip 7.11317pt\forall y\in\mathbb{R}^{e},$

for $i,j\in\\{1,\cdots,d\\}$. Then

(55) $\displaystyle\sum_{i,j=1}^{d}g_{i}^{\hskip
0.7113pt\prime}(y)g_{j}(y)I_{ji}=\frac{1}{2}\sum_{i,j=1}^{d}g_{i}^{\hskip
0.7113pt\prime}(y)g_{j}(y)\big{(}I_{i}\cdot I_{j}\big{)},\hskip 133.72795pt$
(56) $\displaystyle\hskip 1.42262pt\sum_{i,j,k=1}^{d}\big{(}g_{i}^{\hskip
0.7113pt\prime\prime}(y)\big{(}g_{j}(y),g_{k}(y)\big{)}+g_{i}^{\hskip
0.7113pt\prime}(y)g_{j}^{\hskip
0.7113pt\prime}(y)g_{k}(y)\big{)}I_{kji}=\frac{1}{6}\sum_{i,j,k=1}^{d}\big{(}\cdots\big{)}\big{(}I_{i}\cdot
I_{j}\cdot I_{k}\big{)},$ (57) $\displaystyle\sum_{i,j,k,l=1}^{d}\Big{(}\hskip
0.7113ptg_{i}^{\hskip
0.7113pt\prime\prime\prime}(y)\big{(}g_{j}(y),g_{k}(y),g_{l}(y)\big{)}+g_{i}^{\hskip
0.7113pt\prime\prime}(y)\big{(}g_{j}^{\hskip
0.7113pt\prime}(y)g_{l}(y),g_{k}(y)\big{)}\hskip 72.55452pt$
$\displaystyle+\,\,\,g_{i}^{\hskip
0.7113pt\prime\prime}(y)\big{(}g_{j}^{\hskip
0.7113pt\prime}(y)g_{k}(y),g_{l}(y)\big{)}+g_{i}^{\hskip
0.7113pt\prime\prime}(y)\big{(}g_{j}(y),g_{k}^{\hskip
0.7113pt\prime}(y)g_{l}(y)\big{)}\hskip 56.9055pt$
$\displaystyle+\,\,\,g_{i}^{\hskip 0.7113pt\prime}(y)g_{j}^{\hskip
0.7113pt\prime\prime}(y)\big{(}g_{k}(y),g_{l}(y)\big{)}+g_{i}^{\hskip
0.7113pt\prime}(y)g_{j}^{\hskip 0.7113pt\prime}(y)g_{k}^{\hskip
0.7113pt\prime}(y)g_{l}(y)\Big{)}I_{lkji}\hskip 48.36967pt$
$\displaystyle=\frac{1}{24}\sum_{i,j,k=1}^{d}\big{(}\cdots\big{)}\big{(}I_{i}\cdot
I_{j}\cdot I_{k}\cdot I_{l}\big{)},\hskip 106.69783pt$

where $(\hskip 0.7113pt\cdots)$ denote the same sum of vector field
derivatives as in the left hand sides.

######  Proof A.10.

Equation (55) follows directly by the commutativity condition (54) and the
antisymmetry of the $I_{[i,j]}$ term in (51). Differentiating both sides of
(54) gives

$\displaystyle g_{i}^{\hskip
0.7113pt\prime\prime}(y)\big{(}g_{j}(y),x\big{)}+g_{i}^{\hskip
0.7113pt\prime}(y)g_{j}^{\hskip 0.7113pt\prime}(y)(x)=g_{j}^{\hskip
0.7113pt\prime\prime}(y)\big{(}g_{i}(y),x\big{)}+g_{j}^{\hskip
0.7113pt\prime}(y)g_{i}^{\hskip 0.7113pt\prime}(y)(x),$

and

$\displaystyle g_{i}^{\hskip
0.7113pt\prime\prime\prime}(y)\big{(}g_{j}(y),x,z\big{)}+g_{i}^{\hskip
0.7113pt\prime\prime}(y)\big{(}g_{j}^{\prime}(y)z,x\big{)}+g_{i}^{\hskip
0.7113pt\prime\prime}(y)\big{(}g_{j}^{\prime}(y)x,z\big{)}+g_{i}^{\hskip
0.7113pt\prime}(y)g_{j}^{\hskip 0.7113pt\prime\prime}(y)(x,z)$
$\displaystyle\hskip 14.22636pt=g_{j}^{\hskip
0.7113pt\prime\prime\prime}(y)\big{(}g_{i}(y),x,z\big{)}+g_{j}^{\hskip
0.7113pt\prime\prime}(y)\big{(}g_{i}^{\prime}(y)z,x\big{)}+g_{j}^{\hskip
0.7113pt\prime\prime}(y)\big{(}g_{i}^{\prime}(y)x,z\big{)}+g_{j}^{\hskip
0.7113pt\prime}(y)g_{i}^{\hskip 0.7113pt\prime\prime}(y)(x,z),$

for $x,y,z\in\mathbb{R}^{e}$. These identities imply that the $(\hskip
0.7113pt\cdots)$ terms in (56) and (57) are symmetric in the $(i,j)$ indices.
In addition, the terms in (56) are $(j,k)$-symmetric as

$\displaystyle(\,\cdots)=\underbrace{g_{i}^{\hskip
0.7113pt\prime\prime}(y)}_{\begin{subarray}{c}\text{symmetric }\\\ \text{and
bilinear}\end{subarray}}\big{(}g_{j}(y),g_{k}(y)\big{)}+g_{i}^{\hskip
0.7113pt\prime}(y)\underbrace{g_{j}^{\hskip
0.7113pt\prime}(y)g_{k}(y)}_{=\,g_{k}^{\hskip 0.7113pt\prime}(y)g_{j}(y)}.$

By differentiating both sides of $g_{i}^{\hskip
0.7113pt\prime\prime}(g_{j},g_{k})+g_{i}^{\hskip 0.7113pt\prime}g_{j}^{\hskip
0.7113pt\prime}g_{k}=g_{i}^{\hskip
0.7113pt\prime\prime}(g_{k},g_{j})+g_{i}^{\hskip 0.7113pt\prime}g_{k}^{\hskip
0.7113pt\prime}g_{j}\hskip 0.7113pt$, we obtain

$\displaystyle g_{i}^{\hskip
0.7113pt\prime\prime\prime}(y)\big{(}g_{j}(y),g_{k}(y),x\big{)}+g_{i}^{\hskip
0.7113pt\prime\prime}(y)\big{(}g_{j}^{\hskip
0.7113pt\prime}(y)x,g_{k}(y)\big{)}+g_{i}^{\hskip
0.7113pt\prime\prime}(y)\big{(}g_{j}(y),g_{k}^{\prime}(y)x\big{)}$
$\displaystyle\hskip 14.22636pt+g_{i}^{\hskip
0.7113pt\prime\prime}(y)\big{(}g_{j}^{\hskip
0.7113pt\prime}(y)g_{k}(y),x\big{)}+g_{i}^{\hskip
0.7113pt\prime}(y)g_{j}^{\hskip
0.7113pt\prime\prime}(y)\big{(}g_{k}(y),x\big{)}+g_{i}^{\hskip
0.7113pt\prime}(y)g_{j}^{\hskip 0.7113pt\prime}(y)g_{k}^{\hskip
0.7113pt\prime}(y)x$ $\displaystyle\hskip 42.67912pt=g_{i}^{\hskip
0.7113pt\prime\prime\prime}(y)\big{(}g_{k}(y),g_{j}(y),x\big{)}+g_{i}^{\hskip
0.7113pt\prime\prime}(y)\big{(}g_{k}^{\hskip
0.7113pt\prime}(y)x,g_{j}(y)\big{)}+g_{i}^{\hskip
0.7113pt\prime\prime}(y)\big{(}g_{k}(y),g_{j}^{\prime}(y)x\big{)}$
$\displaystyle\hskip 71.13188pt+g_{i}^{\hskip
0.7113pt\prime\prime}(y)\big{(}g_{k}^{\hskip
0.7113pt\prime}(y)g_{j}(y),x\big{)}+g_{i}^{\hskip
0.7113pt\prime}(y)g_{k}^{\hskip
0.7113pt\prime\prime}(y)\big{(}g_{j}(y),x\big{)}+g_{i}^{\hskip
0.7113pt\prime}(y)g_{k}^{\hskip 0.7113pt\prime}(y)g_{j}^{\hskip
0.7113pt\prime}(y)x,$

for $x,y\in\mathbb{R}^{e}$. Thus, it follows that the $(\hskip
0.7113pt\cdots)$ terms in (57) are also $(j,k)$-symmetric. Finally, it is
straightforward to show that the terms in (57) are $(k,l)$-symmetric as,

$\displaystyle(\,\cdots)$ $\displaystyle=\underbrace{g_{i}^{\hskip
0.7113pt\prime\prime\prime}(y)}_{\begin{subarray}{c}\text{symmetric }\\\
\text{and
trilinear}\end{subarray}}\big{(}g_{j}(y),g_{k}(y),g_{l}(y)\big{)}+g_{i}^{\hskip
0.7113pt\prime}(y)\underbrace{g_{j}^{\hskip
0.7113pt\prime\prime}(y)}_{\begin{subarray}{c}\text{symmetric }\\\ \text{and
bilinear}\end{subarray}}\big{(}g_{k}(y),g_{l}(y)\big{)}$ $\displaystyle\hskip
14.22636pt+\underbrace{g_{i}^{\hskip
0.7113pt\prime\prime}(y)\big{(}g_{j}^{\hskip
0.7113pt\prime}(y)g_{l}(y),g_{k}(y)\big{)}+g_{i}^{\hskip
0.7113pt\prime\prime}(y)\big{(}g_{j}^{\hskip
0.7113pt\prime}(y)g_{k}(y),g_{l}(y)\big{)}}_{\text{symmetric in }(k,\hskip
0.7113ptl)}$ $\displaystyle\hskip 14.22636pt+g_{i}^{\hskip
0.7113pt\prime\prime}(y)\big{(}g_{j}(y),\underbrace{g_{k}^{\hskip
0.7113pt\prime}(y)g_{l}(y)}_{=\,g_{l}^{\hskip
0.7113pt\prime}(y)g_{k}(y)}\big{)}+g_{i}^{\hskip
0.7113pt\prime}(y)g_{j}^{\hskip 0.7113pt\prime}(y)\underbrace{g_{k}^{\hskip
0.7113pt\prime}(y)g_{l}(y)}_{=\,g_{l}^{\hskip 0.7113pt\prime}(y)g_{k}(y)}.$

Since the $(\hskip 0.7113pt\cdots)$ terms in (57) have symmetries in the
$(i,j),(j,k)$ and $(k,l)$ indices, they must also be $(i,l)$-symmetric. So by
Theorem A.5, all of the antisymmetric terms in the decompositions of $I_{kji}$
and $I_{lkji}$ will cancel out in their respective sums. Thus, only the
“symmetric parts” of the integrals contribute to the sums (56) and (57). The
result now follows from the decompositions (52) and (53) in Theorem A.5.

##  Appendix B Unbiased approximation of high order iterated integrals

In this section, we derive estimators for certain iterated stochastic
integrals using a polynomial expansion of Brownian motion [29]. We use this
expansion since its first two coefficients give the path’s increment and
space-time Lévy area (Definition 1.1). Just as in [29], the integral that we
would primarily like to approximate is the so-called “space-space-time” Lévy