# An optimal polynomial approximation of
Brownian motion

James Foster††Mathematical Institute, University of Oxford, Woodstock Road,
Oxford, OX2 6GG, UK.
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Research supported by the Engineering and
Physical Sciences Research Council [EP/N509711/1]. Terry Lyons Supported by
the EPSRC grant DATASIG, Alan Turing Institute and Oxford-Man Institute.
Harald Oberhauser††footnotemark:

###### Abstract

In this paper, we will present a strong (or pathwise) approximation of
standard Brownian motion by a class of orthogonal polynomials. The
coefficients that are obtained from the expansion of Brownian motion in this
polynomial basis are independent Gaussian random variables. Therefore it is
practical (requires $N$ independent Gaussian coefficients) to generate an
approximate sample path of Brownian motion that respects integration of
polynomials with degree less than $N$. Moreover, since these orthogonal
polynomials appear naturally as eigenfunctions of the Brownian bridge
covariance function, the proposed approximation is optimal in a certain
weighted $L^{2}(\mathbb{P})$ sense.
In addition, discretizing Brownian paths as piecewise parabolas gives a
locally higher order numerical method for stochastic differential equations
(SDEs) when compared to the piecewise linear approach. We shall demonstrate
these ideas by simulating Inhomogeneous Geometric Brownian Motion (IGBM). This
numerical example will also illustrate the deficiencies of the piecewise
parabola approximation when compared to a new version of the asymptotically
efficient log-ODE (or Castell-Gaines) method.

###### keywords:

Brownian motion, polynomial approximation, numerical methods for SDEs

###### AMS:

41A10, 60J65, 60L90, 65C30

##  1 Introduction

Brownian motion is a central object for modelling real-world systems that
evolve under the influence of random perturbations [1]. In applications where
methods discretize Brownian motion, usually only increments of the path are
generated [2]. In this setting, the best $L^{2}(\mathbb{P})$ approximation of
Brownian motion that is measurable with respect to these increments is given
by the piecewise linear path that agrees on discretization points [3]. This
motivates the following natural question: Are there better discrete
approximations of Brownian motion than piecewise linear? The next simplest
approximant would be a piecewise polynomial, though it is not clear whether
this would be advantageous for tackling problems such as SDE simulation. This
paper can be viewed as a logical continuation of [4], where a polynomial
wavelet representation of Brownian motion was proposed. These wavelets were
constructed to capture certain “geometrical features” of the path, namely the
integrals of the Brownian motion against monomials. We shall investigate the
practical applications of these polynomials and their geometrical features in
the numerical analysis of SDEs.

Fig. 1: Sample paths of Brownian motion with corresponding polynomial
approximations.

The paper is organised as follows. In Section 2, we shall state and prove the
main result of the paper (Theorem 5). This will be a Karhunen-Loève theorem
for the Brownian bridge, where the orthogonal functions used in the
approximation are polynomials. Furthermore, we shall explicitly show that each
basis function is proportional to a shifted $(\alpha,\beta)$-Jacobi polynomial
but with the nonstandard exponents $\alpha=\beta=-1$. This enables us to
construct these orthogonal polynomials using recurrence relations, or as the
difference of two shifted Legendre polynomials whose degrees differ by two.
The resulting polynomial expansion of Brownian motion was independently
discovered by Habermann in [5], where a sharp $L^{2}(\mathbb{P})$ convergence
rate of $O\big{(}\frac{1}{\sqrt{n}}\big{)}$ is established111A Matlab
demonstration can be found at
chebfun.org/examples/stats/RandomPolynomials.html. In Section 3, we shall
investigate some significant consequences of the main theorem.

######  Theorem 1.

Let $W$ denote a standard real-valued Brownian motion on $[\hskip
0.28453pt0,1]$. Let $W^{n}$ be the unique $n$-th degree random polynomial with
a root at $0$ and satisfying

$\displaystyle\int_{0}^{1}u^{k}\,dW_{u}^{n}$
$\displaystyle=\int_{0}^{1}u^{k}\,dW_{u}\hskip 0.35565pt,\hskip
8.53581pt\text{for}\hskip 5.69054ptk=0\hskip 0.7113pt,1\hskip
0.7113pt,\cdots\hskip 0.7113pt,n-1\hskip 0.7113pt.$ (1)

Then $W=\hskip 0.7113ptW^{n}+Z^{n}$ , where $Z^{n}$ is a centered Gaussian
process independent of $W^{n}$.

The above theorem has a simple yet striking conclusion, namely that
polynomials can be unbiased approximants of Brownian motion. In addition, the
first non-trivial case ($n=2$) already has interesting applications within the
numerical analysis of SDEs. One reason is that parabolas can capture the “
space-time area” of Brownian motion.

Fig. 2: Brownian motion can be expressed as a (random) parabola plus
independent noise. Moreover, the approximating parabola has the same increment
and time integral as the original path.

Therefore discretizing Brownian motion using a piecewise parabola gives a
locally high order methodology for numerically solving one-dimensional SDEs.
However, since certain triple iterated integrals of Brownian motion and time
are partially matched by these parabolas, we expect this method to have only
an $O(h)$ rate of convergence (where $h$ denotes the step size used). This
gives motivation for the following theorem:

######  Theorem 2.

Let $\wideparen{W}$ be the (unique) quadratic polynomial with a root at $0$
and

$\displaystyle\wideparen{W}_{1}=W_{1}\hskip 0.7113pt,\hskip
14.22636pt\int_{0}^{1}\wideparen{W}_{u}\,du=\int_{0}^{1}W_{u}\,du\hskip
0.7113pt.$

Then the following third order iterated integral of Brownian motion can be
estimated:

$\displaystyle\mathbb{E}\left[\hskip
0.7113pt\int_{0}^{1}W_{u}^{2}\,du\,\Big{|}\,W_{1}\hskip
0.7113pt,\int_{0}^{1}W_{u}\,du\hskip
0.7113pt\right]=\int_{0}^{1}\wideparen{W}_{u}^{2}\,du+\frac{1}{15}\hskip
0.7113pt.$ (2)

The above theorem can be directly incorporated into the stochastic Taylor
method as well as the log-ODE or Castell-Gaines method (see [6], [7]). We will
show that by estimating this non-trivial iterated integral with its
conditional expectation, we can design numerical methods that enjoy high
orders of both strong and weak convergence. Specifically, for a general SDE
that is driven by a one-dimensional Brownian motion and governed by
sufficiently regular vector fields (smooth with bounded derivatives), the
numerical methods that correctly utilize the above conditional expectation
will have a strong convergence rate of $O(h^{\frac{3}{2}})$ as well as a weak
convergence rate of $O(h^{2})$.

These high orders of convergence can also be achieved in the multidimensional
setting provided the vector fields governing the SDE satisfy certain
commutativity conditions. For example, this estimator has applications for
simulating SDEs with additive noise:

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

where $f$ is a smooth vector field on $\mathbb{R}^{d}$, $\sigma>0$ is constant
and $W$ is now $d$-dimensional. By considering Theorem 2, we expect that
$y_{t}$ is well approximated (for small $t$) by

$\displaystyle\wideparen{y}_{t}+\frac{1}{30}\hskip
0.7113ptt^{2}\sigma^{2}\Delta f(y_{0})\hskip 0.7113pt,$

where $\wideparen{y}$ denotes the solution of the below ODE driven by a “
Brownian parabola” $\wideparen{W}$,

$\displaystyle d\hskip 0.7113pt\wideparen{y}_{t}$ $\displaystyle=f(\hskip
0.7113pt\wideparen{y}_{t})\,dt+\sigma\hskip 0.7113ptd\wideparen{W}_{t}\hskip
0.7113pt,$ $\displaystyle\wideparen{y}_{0}$ $\displaystyle=y_{0}\hskip
0.7113pt.$

This parabola-driven ODE can then be discretized using a three-stage Runge-
Kutta method and the resulting SDE approximation shall be investigated in a
future work.

Since these methods are based on the conditional expectation given by Theorem
2, they are designed to minimize the leading error term within local Taylor
expansions. This sense of optimality is conceptually similar to that of the
asymptotically efficient SDE approximations developed by Clark [8], Newton [9,
10] and Castell & Gaines [6]. The key difference is that we are employing
additional integral information about $W$. Hence, this line of research could
provide a further insight into the approximation of Itô integrals using linear
path information, where there already are a number of results concerning the
computational complexity of methods (see [11], [12] and [13]).
Most notably, Tang and Xiao [11] consider the same triple iterated integral as
in (2) and present an asymptotically optimal approximation that performs well
when a limited number of random variables are used (see Table 2 for these
numerical results). Whilst there are other senses of optimality (such as those
discussed in [14] and [15]) that could be used when analysing the proposed
approximations of Brownian motion and SDE solutions, we shall estimate errors
in an $L^{2}(\mathbb{P})$ sense throughout the paper. In particular, we can
apply the main result to quantify the error of the new estimator.

######  Theorem 3.

Using the same notation as before, we have the following variance:

$\displaystyle\operatorname{Var}\left(\int_{0}^{1}W_{u}^{2}\,du\,\Big{|}\,W_{1}\hskip
0.7113pt,\int_{0}^{1}W_{u}\,du\right)=\frac{11}{6300}+\frac{1}{180}\hskip
0.7113ptW_{1}^{2}+\frac{1}{175}\left(\hskip
0.7113pt\int_{0}^{1}W_{u}\,du-\frac{1}{2}\hskip 0.7113ptW_{1}\right)^{2}.$

In Section 4, we demonstrate the applicability of these ideas to SDE
simulation through various discretizations of Inhomogeneous Geometric Brownian
Motion (IGBM)

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

where $a\geq 0$ and $b\in\mathbb{R}$ are the mean reversion parameters and
$\sigma\geq 0$ is the volatility.

In mathematical finance, IGBM is an example of a short rate model that can be
both mean-reverting and non-negative. It is therefore suitable for modelling
interest rates, stochastic volatilities and default intensities [16]. From a
mathematical viewpoint, IGBM is one of the simplest SDEs that has no known
method of exact simulation [17]. By incorporating the ideas provided by the
main theorem into the log-ODE method, we will produce a state-of-the-art
numerical approximation of IGBM. Although the vector fields for IGBM are not
bounded, our numerical evidence indicates that the method has a strong
convergence rate of $O(h^{\frac{3}{2}})$ and a weak convergence rate of
$O(h^{2})$.

###  1.1 Notation

Below is some of the notation that is used throughout the paper.

Symbol |  Meaning | Page
---|---|---
$W$ |  a standard real-valued Brownian motion. |  2
$B$ |  a standard real-valued Brownian bridge on $[\hskip 0.28453pt0\hskip 0.14226pt,1\hskip 0.28453pt]$. |  5
$\mu$ |  a Borel measure on $[\hskip 0.28453pt0,1]$ defined by a singular weight function. |  5
|  $\begin{aligned} \hskip 85.35826pt\mu(a,b)=\int_{a}^{b}\frac{1}{x(1-x)}\,dx,\end{aligned}$ |
|  for all open intervals $(a,b)\subset[\hskip 0.28453pt0,1]$. |
$\left\\{e_{k}\right\\}_{k\geq 1}$ |  a family of Jacobi-like polynomials with $\deg\left(e_{k}\right)=k+1$ that are orthogonal with respect to weight function $w(x):=\hskip 0.7113pt\frac{1}{x(1-x)}\hskip 0.7113pt.$ |  5
$I_{k}$ |  a time integral of $B$ times the polynomial $e_{k}(t)\,w(t)$ over $[\hskip 0.28453pt0\hskip 0.14226pt,1\hskip 0.28453pt]$, |  5
|  $\begin{aligned} \hskip 85.35826ptI_{k}=\int_{0}^{1}B_{t}\cdot\frac{e_{k}\left(t\right)}{t(1-t)}\,dt.\end{aligned}$ |
$K_{B}$ |  the covariance function of $B$, that is $K_{B}(s,t)=\min(s,t)-st$. |  5
$P_{k}^{(\alpha,\beta)}$ |  the $k$-th order $(\alpha,\beta)$-Jacobi polynomial on $[-1,1\hskip 0.28453pt]$ $\left(\alpha,\beta>-1\right)$. | 11
$Q_{k}$ |  the $k$-th order Legendre polynomial on $[-1,1\hskip 0.28453pt]$, i.e. $Q_{k}=P_{k}^{(0,0)}$. | 13
$y$ |  a solution of the Stratonovich SDE on the finite interval $[\hskip 0.28453pt0\hskip 0.14226pt,T\hskip 0.28453pt]$, | 14
|  $\begin{aligned} \hskip 73.97716ptdy_{t}&=f_{0}(y_{t})\,dt+f_{1}(y_{t})\circ dW_{t}\hskip 0.7113pt,\\\\[3.0pt] y_{0}&=\xi,\end{aligned}$ |
|  where $y,\xi\in\mathbb{R}^{e}$, and $f_{i}:\mathbb{R}^{e}\rightarrow\mathbb{R}^{e}$ denote smooth vector fields. |
|  (Itô SDEs will be defined on fixed intervals with the same form) |
$[s,t]$ |  a general closed subinterval of $[\hskip 0.42677pt0\hskip 0.14226pt,T\hskip 0.28453pt]$, usually considered small. | 14
$h$ |  the step size that a numerical method uses, typically $h=t-s$. | 14
$W_{s,t}$ |  the increment of Brownian motion over $[\hskip 0.28453pts\hskip 0.14226pt,t\hskip 0.28453pt]$, $W_{s,t}:=W_{t}-W_{s}\hskip 0.14226pt$. | 14
$\wideparen{W}$ |  the Brownian parabola corresponding to $W$ over some interval. | 15
|  $\begin{aligned} \hskip 9.95845pt\wideparen{W}_{u}=W_{s}+\frac{u-s}{h}\,W_{s,t}+\frac{6\hskip 0.7113pt(u-s)(t-u)}{h^{2}}\,H_{s,t}\hskip 0.7113pt,\hskip 8.25128pt\forall u\in[\hskip 0.28453pts\hskip 0.14226pt,t\hskip 0.28453pt].\end{aligned}$ |
$Z$ |  the Brownian arch corresponding to $W$ defined as $Z:=W-\wideparen{W}$. | 15
$H_{s,t}$ |  the rescaled space-time Lévy area of Brownian motion on $[\hskip 0.28453pts\hskip 0.14226pt,t\hskip 0.28453pt]$, | 15
|  $\begin{aligned} \hskip 64.01869ptH_{s,t}&=\frac{1}{h}\int_{s}^{t}W_{s,u}-\frac{u-s}{h}\,W_{s,t}\,du.\end{aligned}$ |
$L_{s,t}$ |  the space-space-time Lévy area of Brownian motion over $[\hskip 0.28453pts\hskip 0.14226pt,t\hskip 0.28453pt]$, | 17
|  $\begin{aligned} L_{s,t}&=\frac{1}{6}\left(\int_{s}^{t}\int_{s}^{u}\int_{s}^{v}\circ\,dW_{r}\circ dW_{v}\,du-2\int_{s}^{t}\int_{s}^{u}\int_{s}^{v}\circ\,dW_{r}\,dv\circ dW_{u}\right.\\\ &\hskip 28.45274pt\left.+\int_{s}^{t}\int_{s}^{u}\int_{s}^{v}\,dr\circ dW_{v}\circ dW_{u}\right),\end{aligned}$ |
$Y$ |  an approximation for the true solution $y$ of a Stratonovich SDE. | 19
$[\,\cdot\hskip 1.42262pt,\cdot\,]$ |  the standard Lie bracket of vector fields, $[f_{0},f_{1}]=f_{1}^{\prime}\hskip 0.35565ptf_{0}-f_{0}^{\prime}\hskip 0.35565ptf_{1}$. | 19

##  2 Main result

It was shown in [4] that Brownian motion can be generated using Alpert-Rokhlin
multiwavelets (see [18]). The mother functions that generate this wavelet
basis are supported on $[\hskip 0.28453pt0,1]$ and are defined using
polynomials as follows:

######  Definition 4 (Alpert-Rokhlin wavelets).

For $q\geq 1$, define the $q$ functions $\phi^{q,1},\cdots,\phi^{q,q}:[\hskip
0.42677pt0,1\hskip 0.42677pt]\rightarrow\mathbb{R}$ as piecewise polynomials
of degree $q-1$ with pieces on $[\hskip 0.42677pt0,\frac{1}{2}\hskip
0.42677pt]$, $[\hskip 0.42677pt\frac{1}{2},1\hskip 0.42677pt]$ that satisfy
the following conditions for all $p\in\left\\{1,\,\cdots\,,q\hskip
0.42677pt\right\\}$ and $t\in[\hskip 0.42677pt0,\frac{1}{2}\hskip 0.42677pt)$
$:$

$\displaystyle\phi^{q,p}(t)$ $\displaystyle=(-1)^{q+p-1}\phi^{q,p}(1-t),$ (3)
$\displaystyle\int_{0}^{1}\phi^{q,p}(t)\phi^{q,r}(t)\,dt$
$\displaystyle=\delta_{qr}\hskip 0.7113pt,\hskip 8.53581pt\text{for}\hskip
5.69054pt1\leq r\leq q,$ (4) $\displaystyle\int_{0}^{1}t^{k}\phi^{q,p}(t)\,dt$
$\displaystyle=0,\hskip 16.07578pt\text{for}\hskip 6.25958pt0\leq k\leq q-1.$
(5)

The Alpert-Rokhlin multiwavelets of order $q$ can now be generated by
translating and scaling the mother functions $\phi^{q,p}$.

$\displaystyle\phi_{nk}^{q,p}(t):=\frac{1}{\sqrt{2^{n}}}\,\phi^{q,p}(2^{n}t-k),$

for $n\geq 0$ and $k\in\left\\{0,\cdots,2^{n}-1\right\\}$.

Whilst our results will not be presented in terms of the above wavelets, we
shall see that the polynomials of interest are directly related to conditions
(3), (4) and (5). The main result of this paper gives an effective method for
approximating sample paths of Brownian motion by a class of Jacobi-like
polynomials. The proof is based on the interpretation of these polynomials as
eigenfunctions of an integral operator defined by the Brownian bridge
covariance function222The Brownian bridge is the centered Gaussian process
with covariance $K_{B}(s,t)=\min(s,t)-st$.. These orthogonal polynomials,
which lie at the heart of this paper, will also help us interpret the
geometrical features that certain normally distributed iterated integrals
encode about the Brownian path.

######  Theorem 5 (A polynomial Karhunen-Loève theorem for the Brownian
bridge).

Let $B$ denote a Brownian bridge on $[\hskip 0.28453pt0,1]$ and consider the
Borel measure $\mu$ given by

$\displaystyle\mu(a,b):=\int_{a}^{b}\frac{1}{x(1-x)}\,dx,\hskip
9.95845pt\text{for all open intervals}\hskip 2.84526pt(a,b)\subset[\hskip
0.28453pt0,1].$

Then there exists a family of orthogonal polynomials $\\{e_{k}\\}_{k\geq 1}$
with $\deg\left(e_{k}\right)=k+1$ and

$\displaystyle\int_{0}^{1}e_{i}\,e_{j}\,d\mu=\delta_{ij}\hskip 0.35565pt,$

with $\delta_{ij}$ denoting the Kronecker delta, such that $B$ admits the
following representation

$\displaystyle B=\sum_{k=1}^{\infty}I_{k}e_{k}\hskip 0.35565pt,$ (6)

where $\\{I_{k}\\}$ is the collection of independent centered Gaussian random
variables with

$\displaystyle
I_{k}:=\int_{0}^{1}B_{t}\cdot\frac{e_{k}\left(t\right)}{t(1-t)}\,dt,$ (7)

and

$\displaystyle\operatorname{Var}(I_{k})=\frac{1}{k(k+1)}\,.$

Furthermore, $\left\\{e_{k}\right\\}$ is an optimal orthonormal basis of
$L^{2}([\hskip 0.28453pt0,1],\mu)$ for approximating $B$ by truncated series
expansions with respect to the following weighted $L^{2}(\mathbb{P})$ norm

$\displaystyle\left\|X\right\|_{L_{\mu}^{2}(\mathbb{P})}:=\sqrt{\hskip
0.7113pt\mathbb{E}\left[\hskip
0.7113pt\int_{0}^{1}\left(X_{s}\right)^{2}\,d\mu(s)\right]},$

where $X$ is a square $\mu$-integrable process.

###### Proof.

Our argument is that of the Karhunen-Loève theorem in general $L^{2}$ spaces.
Note that $B$ is a square $\mu$-integrable process as

$\displaystyle\mathbb{E}\left[\hskip
0.7113pt\int_{0}^{1}(B_{s})^{2}\,d\mu(s)\right]=\int_{0}^{1}\mathbb{E}\big{[}(B_{s})^{2}\big{]}\,d\mu(s)=\int_{0}^{1}s(1-s)\cdot\frac{1}{s(1-s)}\,ds=1<\infty.$

Let $K_{B}$ denote the covariance function for the standard Brownian bridge on
$[\hskip 0.28453pt0,1]$. Since $K_{B}(s,t)=\min(s,t)-st$, it can be shown by
direct calculation that $K_{B}$ satisfies

$\displaystyle\left\|K_{B}\right\|_{L^{2}([\hskip
0.28453pt0,1]^{2},\,\mu^{2})}^{2}=\int_{0}^{1}\int_{0}^{1}(\min(s,t)-st)^{2}\,d\mu(s)d\mu(t)=\frac{1}{3}\pi^{2}-3<\infty.$

Hence, it follows that the integral operator $T_{K}:L^{2}([\hskip
0.28453pt0,1],\,\mu)\rightarrow L^{2}([\hskip 0.28453pt0,1],\,\mu)$ given by

$\displaystyle(T_{K}f)(t):=\int_{0}^{1}K_{B}(s,t)f(s)\,d\mu(s),$

is well-defined and continuous. In addition, the variance function
$k_{B}(x):=K_{B}(x,x)$ for $x\in[\hskip 0.28453pt0,1]$ is $\mu$-integrable as

$\displaystyle\int_{0}^{1}|k_{B}(x)|\,d\mu(x)=\int_{0}^{1}x(1-x)\cdot\frac{1}{x(1-x)}\,dx=1<\infty.$

Therefore, we can apply Mercer’s theorem for kernels on general $L^{2}$ spaces
(see [19]). It then follows from Mercer’s theorem that there exists an
orthonormal set $\\{e_{k}\\}_{k\geq 1}$ of $L^{2}([\hskip
0.28453pt0,1],\,\mu)$ consisting of eigenfunctions of $T_{K}$ such that the
corresponding sequence of eigenvalues $\\{\lambda_{k}\\}_{k\geq 1}$ is non-
negative. Moreover, the eigenfunctions corresponding to non-zero eigenvalues
are continuous on $[\hskip 0.28453pt0,1]$ and the kernel $K_{B}$ has the
representation

$\displaystyle K_{B}(s,t)=\sum_{k=1}^{\infty}\lambda_{k}e_{k}(s)e_{k}(t),$ (8)

where the series (8) converges absolutely and uniformly on compact subsets of
$[\hskip 0.28453pt0,1]$.

In the next part of the proof, we will see that each $e_{k}$ is a polynomial
of degree $k+1$. As each $e_{k}$ is an eigenfunction of $T_{K}$, we have

$\displaystyle\int_{0}^{1}\frac{\min(s,t)-st}{s(1-s)}\,e_{k}(s)\,ds=\lambda_{k}e_{k}(t).$
(9)

Since $e_{k}\in L^{2}([\hskip 0.28453pt0,1],\,\mu)$, it follows that
$e_{k}(0)=0$ and $e_{k}(1)=0$ for each $k\geq 1$. Therefore by using the
Leibniz integral rule to twice differentiate both sides of (9) and then
multiplying by $t(1-t)$, we observe that $e_{k}$ satisfies the differential
equation

$\displaystyle t(1-t)\lambda_{k}e_{k}^{\prime\prime}(t)+e_{k}(t)=0.$ (10)

Since $e_{k}\neq 0$, we have that $\lambda_{k}\neq 0$. Differentiating the LHS
of the ODE (10) produces

$\displaystyle
t(1-t)\frac{d^{2}}{dt^{2}}(e_{k}^{\prime})+(1-2t)\frac{d}{dt}(e_{k}^{\prime})+\frac{1}{\lambda_{k}}e_{k}^{\prime}=0.$

For $x\in[-1,1]$, we define the function

$\displaystyle y_{k}(x):=e_{k}^{\prime}\left(\frac{1}{2}(1+x)\right).$

Thus $y_{k}$ satisfies the following differential equation

$\displaystyle(1-x^{2})y_{k}^{\prime\prime}(x)-2xy_{k}^{\prime}(x)+\frac{1}{\lambda_{k}}y_{k}(x)=0.$
(11)

Remarkably, this is the Legendre differential equation [20]. It then follows
using classical Sturm-Liouville theory that $\frac{1}{\lambda_{k}}=k(k+1)$ and
$y_{k}$ is proportional to the $k$-th Legendre polynomial. Therefore, the
derivative $e_{k}^{\prime}$ will be a constant multiple of the $k$-th shifted
Legendre polynomial and hence each $e_{k}$ is a polynomial of degree $k+1$.

We can now define the following integrals for $k\geq 1$,

$\displaystyle
I_{k}:=\int_{0}^{1}B_{t}\cdot\frac{e_{k}\left(t\right)}{t(1-t)}\,dt.$

It follows from Fubini’s theorem that

$\displaystyle\mathbb{E}[I_{k}]$ $\displaystyle=0,$
$\displaystyle\mathbb{E}[I_{i}I_{j}]$ $\displaystyle=\mathbb{E}\left[\hskip
0.7113pt\int_{0}^{1}\int_{0}^{1}B_{s}B_{t}\,e_{i}(s)\,e_{j}(t)\,d\mu(s)\,d\mu(t)\right]$
$\displaystyle=\int_{0}^{1}\int_{0}^{1}\mathbb{E}[B_{s}B_{t}]\,e_{i}(s)\,e_{j}(t)\,d\mu(s)\,d\mu(t)$
$\displaystyle=\int_{0}^{1}\,e_{j}(t)\left(\int_{0}^{1}K_{B}(s,t)\,e_{i}(s)\,d\mu(s)\right)\,d\mu(t)$
$\displaystyle=\lambda_{i}\delta_{ij}.$

Since each $I_{k}$ is defined by a linear functional on the same Gaussian
process $B$, we see from the above that $\\{I_{k}\\}$ is a collection of
uncorrelated (and therefore independent) Gaussian random variables with

$\displaystyle\mathbb{E}[I_{k}]$ $\displaystyle=0,$
$\displaystyle\operatorname{Var}(I_{k})$ $\displaystyle=\frac{1}{k(k+1)}.$

Finally, the $L^{2}(\mathbb{P})$ convergence we require follows as

$\displaystyle\mathbb{E}\left[\left(B_{t}-\sum_{k=1}^{N}I_{k}e_{k}(t)\right)^{2}\right]$
$\displaystyle=k_{B}(t)+\mathbb{E}\left[\hskip
0.7113pt\sum_{i,j=1}^{N}I_{i}I_{j}\,e_{i}(t)\hskip
0.7113pte_{j}(t)\right]-2\hskip 0.7113pt\mathbb{E}\left[\hskip
0.7113ptB_{t}\sum_{k=1}^{N}I_{k}e_{k}(t)\right]$
$\displaystyle=k_{B}(t)+\sum_{i=1}^{N}\lambda_{i}e_{i}^{2}(t)-2\hskip
0.7113pt\mathbb{E}\left[\hskip
0.7113pt\sum_{i=1}^{N}\int_{0}^{1}B_{s}B_{t}\,e_{i}(s)\hskip
0.7113pte_{i}(t)\,d\mu(s)\right]$
$\displaystyle=k_{B}(t)-\sum_{i=1}^{N}\lambda_{i}e_{i}^{2}(t),$

which converges to $0$ by Mercer’s theorem (8).

All that remains is to prove optimality for the truncated series expansions of
(6). Let $\left\\{f_{k}\right\\}_{k\geq 1}$ denote an orthonormal basis of
$L^{2}([\hskip 0.28453pt0,1],\mu)$ such that

$\displaystyle B=\sum_{k=1}^{\infty}J_{k}f_{k},\hskip
8.53581pt\text{where}\hskip
8.53581ptJ_{k}:=\int_{0}^{1}B_{t}\,f_{k}(t)\,d\mu(t),\hskip 8.53581pt\forall
k\geq 1.$

For $n\geq 1$, we consider an error process associated with the above:
$r_{n}:=\sum\limits_{k=n+1}^{\infty}J_{k}f_{k}\,$.

Then the square $L^{2}(\mathbb{P})$ norm of the $n$-th error process admits
the following expansion,

$\displaystyle\left\|r_{n}(t)\right\|_{L^{2}(\mathbb{P})}^{2}$
$\displaystyle=\mathbb{E}\left[\hskip
0.7113pt\sum_{i=n+1}^{\infty}\sum_{j=n+1}^{\infty}J_{i}J_{j}f_{i}(t)f_{j}(t)\right]$
$\displaystyle=\sum_{i=n+1}^{\infty}\sum_{j=n+1}^{\infty}\mathbb{E}\left[\hskip
0.7113pt\int_{0}^{1}\int_{0}^{1}B_{s}B_{t}\,f_{i}(s)f_{j}(t)\,d\mu(s)\,d\mu(t)\right]f_{i}(t)f_{j}(t)$
$\displaystyle=\sum_{i=n+1}^{\infty}\sum_{j=n+1}^{\infty}\left(\int_{0}^{1}\int_{0}^{1}K_{B}(s,t)\,f_{i}(s)f_{j}(t)\,d\mu(s)\,d\mu(t)\right)f_{i}(t)f_{j}(t).$

Integrating the above with respect to $\mu$ and using the orthogonality of
$\\{f_{k}\\}_{k\geq 1}$ gives

$\displaystyle\left\|r_{n}\right\|_{L_{\mu}^{2}(\mathbb{P})}^{2}=\int_{0}^{1}\left\|r(t)\right\|_{L^{2}(\mathbb{P})}^{2}\,d\mu(t)=\sum_{k=n+1}^{\infty}\int_{0}^{1}\int_{0}^{1}K_{B}(s,t)f_{k}(s)f_{k}(t)\,d\mu(s)\,d\mu(t).$

Note that any optimal orthonormal basis of $L^{2}([\hskip 0.28453pt0,1],\mu)$
solves the following problem:

$\displaystyle\min_{f_{k}}\left\|r_{n}\right\|_{L_{\mu}^{2}(\mathbb{P})}^{2}\hskip
8.53581pt\text{subject to}\hskip 8.53581pt\left\|f_{k}\right\|_{L^{2}([\hskip
0.28453pt0,1],\mu)}=1.$

By introducing Lagrange multipliers $\nu_{k}$, we wish to find functions
$\\{f_{k}\\}$ that minimize

$\displaystyle
E_{n}[\left\\{f_{k}\right\\}]:=\sum_{k=n+1}^{\infty}\int_{0}^{1}\int_{0}^{1}K_{B}(s,t)f_{k}(s)f_{k}(t)d\mu(s)d\mu(t)-\nu_{k}\left(\int_{0}^{1}\big{(}f_{k}(s)\big{)}^{2}d\mu(s)-1\right).$

We will now consider the following square integrable functions, defined for
$s,t\in\left(0,1\right)$:

$\displaystyle\tilde{f}_{k}(t):=f_{k}(t)\cdot\frac{1}{\sqrt{t(1-t)}}\,,\hskip
14.22636pt\tilde{K}_{B}(s,t):=K_{B}(s,t)\cdot\frac{1}{\sqrt{s(1-s)}}\cdot\frac{1}{\sqrt{t(1-t)}}\,.$

Therefore it is enough to find a family of functions $\\{\tilde{f}_{k}\\}$ in
$L^{2}([\hskip 0.28453pt0,1])$ which minimizes

$\displaystyle\tilde{E}_{n}[\\{\tilde{f}_{k}\\}]:=\sum_{k=n+1}^{\infty}\int_{0}^{1}\int_{0}^{1}\tilde{K}_{B}(s,t)\tilde{f}_{k}(s)\tilde{f}_{k}(t)\,ds\,dt-\nu_{k}\left(\int_{0}^{1}\big{(}\tilde{f}_{k}(s)\big{)}^{2}\,ds-1\right).$

To find a minimizer, we set the functional derivative of $\tilde{E}_{n}$ with
respect to $\tilde{f}_{k}$ to zero.

$\displaystyle\frac{\partial\tilde{E}_{n}}{\partial\tilde{f}_{k}(t)}=2\int_{0}^{1}\tilde{K}_{B}(s,t)\tilde{f}_{k}(s)\,ds-2\nu_{k}\tilde{f}_{k}(t)=0.$

By using the definitions of $\tilde{f}_{k}$ and $\tilde{K}_{B}$, it is trivial
to show the above is equivalent to

$\displaystyle\int_{0}^{1}K_{B}(s,t)f_{k}(s)\,d\mu(s)=\nu_{k}f_{k}(t),$

which is satisfied if and only if $f_{k}$ are eigenfunctions of $T_{K}$. ∎

This result can naturally be extended to express Brownian motion using
polynomials.

######  Theorem 6.

If $W$ is a standard Brownian motion and $B$ is the associated bridge process
on $[\hskip 0.28453pt0,1]$, then by Theorem 5, we have the below
representation of $W$:

$\displaystyle W=W_{1}e_{0}+\sum_{k=1}^{\infty}I_{k}e_{k}\hskip 0.35565pt,$
(12)

where $e_{0}(t):=t\hskip 0.7113pt$ for $\hskip 0.7113ptt\in[\hskip
0.28453pt0,1]$, and the random variables $\\{I_{k}\\}$ are independent of
$W_{1}$.

Fig. 3: Brownian motion can be expressed as a sum of polynomials with
independent weights. Moreover, these polynomials are orthogonal and capture
different time integrals of the original path.

In the rest of this section, we shall study the key objects introduced in
Theorem 5. Since each orthogonal polynomial lies in $L^{2}([\hskip
0.28453pt0,1],\,\mu)$, it must have roots at $0$ and $1$. Therefore
$e_{k}\cdot\frac{1}{t(1-t)}$ is itself a polynomial but with degree $k-1$, and
one can repeatedly apply the integration by parts formula to the stochastic
integrals $\\{I_{k}\\}$ defined by (7). This enables us to express each
$I_{k}$ in terms of iterated integrals of Brownian motion. Moreover, as
$e_{k}\cdot\frac{1}{t(1-t)}$ has precisely $k-2$ non-zero derivatives, the
highest order iterated integral that is required to fully describe $I_{k}$ is
$\int_{0<s_{1}<\cdots<s_{k}<1}B_{s_{1}}\,ds_{1}\,\cdots\,ds_{k}\hskip
0.35565pt$.

So by applying the integration by parts formula as above, we can construct a
lower triangular $n\times n$ matrix $M_{n}$ with non-zero diagonal entries
that characterizes the relationship between $\\{I_{k}\\}_{1\leq k\leq n}$ and
a set of $n$ iterated integrals of Brownian motion.

Hence, for $n\geq 1$, we can express the $n$ independent Gaussian integrals
$\\{I_{k}\\}_{1\leq k\leq n}$ as

$\displaystyle\begin{pmatrix}I_{1}\\\ \vdots\\\
I_{n}\end{pmatrix}=M_{n}\begin{pmatrix}\int_{0<s_{1}<1}B_{s_{1}}\,ds_{1}\\\
\vdots\\\
\,\int_{0<s_{1}<\cdots<s_{n}<1}B_{s_{1}}\,ds_{1}\,\cdots\,ds_{n}\,\end{pmatrix}.$
(13)

Since $M_{n}$ is an invertible matrix, it follows that the column vectors
appearing in (13) both encode the same information about the Brownian bridge.
This enables us to establish a connection between Brownian motion, iterated
integrals and polynomials.

######  Theorem 7.

Consider the below conditional expectation of Brownian motion,

$\displaystyle W_{t}^{n}:=\mathbb{E}\left[W_{t}\,\Big{|}\,W_{1}\hskip
0.7113pt,\int_{0<s_{1}<1}W_{s_{1}}\,ds_{1}\hskip
0.7113pt,\cdots,\int_{0<s_{1}<\cdots<s_{n-1}<1}W_{s_{1}}\,ds_{1}\cdots\,ds_{n-1}\right].$
(14)

where $t\in[\hskip 0.28453pt0,1]$. Then $W^{n}$ is the unique polynomial of
degree $n$ with a root at $0$ that matches the increment $W_{1}$ and $n-1$
iterated time integrals of the path $\,W$ given by:

$\displaystyle\int_{0<s_{1}<1}W_{s_{1}}\,ds_{1}\hskip
0.7113pt,\cdots,\int_{0<s_{1}<\cdots<s_{n-1}<1}W_{s_{1}}\,ds_{1}\cdots\,ds_{n-1}\hskip
0.35565pt.$ (15)

###### Proof.

It is a direct consequence of (13) that
$W_{t}^{n}=\mathbb{E}[W_{t}\,|\,W_{1},I_{1},\cdots,I_{n-1}]$. Hence by (12)
and independence of the random variables $\\{W_{1},I_{1},\cdots\\}$, we have
that

$\displaystyle W^{n}=W_{1}e_{0}+\sum_{k=1}^{n-1}I_{k}e_{k}\hskip 0.35565pt.$
(16)

Thus $W^{n}$ is indeed a polynomial of degree $n$ with a root at $0$ and that
matches the increment of the Brownian path. Without loss of generality we can
now assume $n\geq 2$. All that remains is to argue $W^{n}$ matches the $n-1$
iterated integrals given in (15). Using the orthogonality of $\\{e_{k}\\}$, it
follows directly from (16) that for $1\leq k\leq n-1$:

$\displaystyle I_{k}$
$\displaystyle=\int_{0}^{1}\Big{(}W_{t}-W_{1}e_{0}(t)\Big{)}\cdot\frac{e_{k}\left(t\right)}{t(1-t)}\,dt$
$\displaystyle=\int_{0}^{1}\left(W_{t}^{n}+\sum_{m=n}^{\infty}I_{m}e_{m}(t)-W_{1}e_{0}(t)\right)\cdot\frac{e_{k}\left(t\right)}{t(1-t)}\,dt$
$\displaystyle=\int_{0}^{1}\Big{(}W_{t}^{n}-W_{1}e_{0}(t)\Big{)}\cdot\frac{e_{k}\left(t\right)}{t(1-t)}\,dt+\underbrace{\sum_{m=n}^{\infty}I_{m}\int_{0}^{1}\frac{e_{m}\left(t\right)e_{k}\left(t\right)}{t(1-t)}\,dt}_{=\,0.}.$

Hence $W^{n}$ matches the integrals of Brownian motion against polynomials
with degree at most $n-1$. By the same argument used in the derivation of
(13), it follows that $W^{n}$ matches the various iterated time integrals
given in the statement of the theorem. The uniqueness of $W^{n}$ is now a
consequence of having $n+1$ different constraints. ∎

###  2.1 Properties of orthogonal polynomials

Although Theorem 5 and Theorem 6 are interesting results from a theoretical
point of view, both lack an explicit construction of the polynomials
$\left\\{e_{k}\right\\}$ that could be implemented in practice. On the other
hand, it was shown that the defining eigenfunction property of each $e_{k}$
implies that its derivative $e_{k}^{\prime}$ is proportional to the $k$-th
shifted Legendre polynomial. Hence the family $\left\\{e_{k}\right\\}$ is the
(normalized) shifted $(\alpha,\beta)$-Jacobi polynomials but with
$\alpha=\beta=-1$. Since Jacobi polynomials are typically studied with
$\alpha,\beta>-1$, it is necessary to show there exists a well-defined limit
when the parameters approach $-1$.

Fig. 4: Sample paths of Brownian motion with corresponding polynomial
approximations.

######  Definition 8.

For $k\geq 2$, the $k$-th degree $(\text{-}1,\text{-}1)$-Jacobi polynomial
$P_{k}^{(\text{-}1,\text{-}1)}$ is

$\displaystyle
P_{k}^{(\text{-}1,\text{-}1)}:=\lim_{\alpha,\beta\rightarrow-1^{+}}P_{k}^{(\alpha,\beta)}.$

Naturally, for this definition to be unambiguous, we will require the
following lemma.

######  Lemma 9.

Let $P_{k}^{(\alpha,\beta)}$ denote the $k$-th degree $(\alpha,\beta)$-Jacobi
polynomial on $[-1,1]$. Then for $k\geq 2$, there exists a real-valued
polynomial $P_{k}$ such that
$\big{\|}P_{k}-P_{k}^{(\alpha,\beta)}\big{\|}_{\infty}\rightarrow 0$ as
$\alpha,\beta\rightarrow-1^{+}$.

###### Proof.

Below is an identity for Jacobi polynomials, with $\alpha,\beta>-1$, given in
[21].

$\displaystyle
P_{k}^{(\alpha,\beta)}(x)=\frac{k+\alpha+\beta+1}{2}\int_{-1}^{x}P_{k-1}^{(\alpha+1,\beta+1)}(u)\,du,\hskip
14.22636pt\text{for all}\hskip 5.69054ptk\geq 2.$ (17)

Therefore, we shall define the $k$-th degree polynomial $P_{k}$ over the
interval $[-1,1]$ by

$\displaystyle
P_{k}(x):=\frac{k-1}{2}\int_{-1}^{x}P_{k-1}^{(0,0)}(u)\,du,\hskip
14.22636pt\text{for all}\hskip 5.69054ptk\geq 2.$ (18)

It is straightforward to verify that $\lim_{\hskip
0.7113pt\alpha,\beta\,\rightarrow\,0}\|\hskip
1.42262ptP_{n}^{(\alpha,\beta)}-\hskip 0.7113ptP_{n}^{(0,0)}\hskip
0.7113pt\|_{\infty}=0$ for $n\in\\{0,1\\}$. So by induction and the recurrence
relation for Jacobi polynomials (see [21]), we have:

$\displaystyle\big{\|}\hskip 1.42262ptP_{n}^{(\alpha,\beta)}-\hskip
0.7113ptP_{n}^{(0,0)}\hskip 0.7113pt\big{\|}_{\infty}\rightarrow 0\hskip
7.11317pt\text{as}\,\,\,\alpha,\beta\,\rightarrow\,0,$

for all $n\geq 0$. Hence by the dominated convergence theorem with (17) and
(18), it follows that $P_{k}^{(\alpha,\beta)}$ will converge pointwise to
$P_{k}$ as $\alpha,\beta\,\rightarrow\,-1$ for each $k\geq 1$. Finally, the
result follows as $P_{k}^{(\alpha,\beta)}$ and $P_{k}$ are always polynomials
with degree $k$. ∎

Using the above definition for $(\text{-}1,\text{-}1)$-Jacobi polynomials, we
can give an explicit formula for the orthonormal polynomials
$\\{e_{k}\\}_{k\geq 1}$ appearing in Theorems 5 and 6.

######  Theorem 10.

Suppose each $e_{k}$ has a positive leading coefficient. Then for $k\geq 1$,

$\displaystyle e_{k}(t)$
$\displaystyle=\frac{1}{k}\sqrt{k(k+1)(2k+1)}\,P_{k+1}^{(\text{-}1,\text{-}1)}(2t-1),\hskip
11.38109pt\forall t\in[\hskip 0.28453pt0,1].$ (19)

###### Proof.

The following identity for $(\alpha,\beta)$-Jacobi polynomials is stated in
[21]:

$\displaystyle\int_{-1}^{1}(1-x)^{\alpha}(1+x)^{\beta}\big{(}P_{n}^{(\alpha,\beta)}(x)\big{)}^{2}\,dx=\frac{2^{\alpha+\beta+1}}{2n+\alpha+\beta+1}\frac{\Gamma(n+\alpha+1)\,\Gamma(n+\beta+1)}{\Gamma(n+\alpha+\beta+1)\,n!}\,,$

for $n\geq 1$ and $\alpha,\beta>-1$. Applying the change of variables,
$t:=\frac{1}{2}(x+1)$, we have

$\displaystyle\int_{0}^{1}t^{\beta}(1-t)^{\alpha}\big{(}P_{n}^{(\alpha,\beta)}(2t-1)\big{)}^{2}dt=\frac{1}{2n+\alpha+\beta+1}\frac{\Gamma(n+\alpha+1)\,\Gamma(n+\beta+1)}{\Gamma(n+\alpha+\beta+1)\,n!}\,,$

for $n\geq 1$ and $\alpha,\beta>-1$. By definition 8, taking the limit
$\alpha,\beta\rightarrow-1^{+}$ will yield

$\displaystyle\int_{0}^{1}\frac{1}{t(1-t)}\left(P_{n}^{(\text{-}1,\text{-}1)}(2t-1)\right)^{2}dt$
$\displaystyle=\frac{1}{2n-1}\frac{(n-1)!\,(n-1)!}{n!\,(n-2)!}$
$\displaystyle=\frac{1}{2n-1}\frac{n-1}{n}\,,\hskip 14.22636pt\text{for
all}\hskip 5.69054ptn\geq 2.$

Therefore by setting $k:=n-1$, we have

$\displaystyle\int_{0}^{1}\frac{1}{t(1-t)}\left(\frac{1}{k}\sqrt{k(k+1)(2k+1)}\,P_{n}^{(\text{-}1,\text{-}1)}(2t-1)\right)^{2}\,dt=1,\hskip
14.22636pt\text{for all}\hskip 3.55658ptk\geq 1.$

Recall that $e_{k}^{\prime}$ is proportional to the $k$-th shifted Legendre
polynomial $P_{k}^{(0,0)}(2t-1)$. Similarly, we saw in the proof of Lemma 9
that the derivative of $P_{k+1}^{(\text{-}1,\text{-}1)}$ is $\frac{k}{2}\hskip
0.7113ptP_{k}^{(0,0)}$. As $e_{k}$ and $P_{k+1}^{(\text{-}1,\text{-}1)}$ are
zero at their respective endpoints, we have that each $e_{k}$ must be
proportional to $P_{k+1}^{(\text{-}1,\text{-}1)}(2t-1)$. The result now
follows from the above calculations. ∎

Having identified an explicit formula for the eigenfunctions $\\{e_{k}\\}$ in
(19), we shall now describe two methodologies for computing the Jacobi-like
polynomials $\big{\\{}P_{k}^{(\text{-}1,\text{-}1)}\big{\\}}$.

The first approach is to use the three-term recurrence relation in the theorem
below.

######  Theorem 11 (Recurrence relation for $(\text{-}1,\text{-}1)$-Jacobi
polynomials).

For $n\geq 2$ ,

$\displaystyle n(n+2)P_{n+2}^{(\text{-}1,\text{-}1)}(x)$
$\displaystyle=(n+1)(2n+1)\hskip
0.7113ptxP_{n+1}^{(\text{-}1,\text{-}1)}(x)-n(n+1)P_{n}^{(\text{-}1,\text{-}1)}(x)\hskip
0.7113pt,$ (20)

where the initial polynomials are given by

$\displaystyle P_{2}^{(\text{-}1,\text{-}1)}(x)$
$\displaystyle=\frac{1}{4}\hskip 0.7113pt(x-1)(x+1)\hskip 0.7113pt,$
$\displaystyle P_{3}^{(\text{-}1,\text{-}1)}(x)$
$\displaystyle=\frac{1}{2}\hskip 0.7113ptx(x-1)(x+1)\hskip 0.7113pt.$

###### Proof.

The below recurrence relation for Jacobi polynomials is presented in [21],

$\displaystyle
2(k+1)(k+\alpha+\beta+1)(2k+\alpha+\beta)P_{k+1}^{(\alpha,\beta)}(x)$
$\displaystyle\hskip
14.22636pt=(2k+\alpha+\beta+1)\left((2k+\alpha+\beta)(2k+\alpha+\beta+2)x+\alpha^{2}-\beta^{2}\right)P_{k}^{(\alpha,\beta)}(x)$
$\displaystyle\hskip
24.18483pt-2(k+\alpha)(k+\beta)(2k+\alpha+\beta+2)\,P_{k-1}^{(\alpha,\beta)}(x)\hskip
0.7113pt,$

for $k\geq 1$ and $\alpha,\beta>-1$. By definition 8, it is possible to take
the limit $\alpha,\beta\rightarrow 1^{+}$ provided that $k\geq 3$. Therefore,
taking this limit and setting $n=k-1\geq 2$ produces

$\displaystyle
4n^{2}(n+2)P_{n+2}^{(\text{-}1,\text{-}1)}(x)=4n(n+1)(2n+1)xP_{n+1}^{(\text{-}1,\text{-}1)}(x)-4n^{2}(n+1)\,P_{n}^{(\text{-}1,\text{-}1)}(x)\hskip
0.7113pt,$

for $n\geq 2$. Dividing the above by $4n$ gives the required recurrence
relation (20). Finally the below formula, stated in [21], can be used to
compute $P_{2}^{(\text{-}1,\text{-}1)}$ and $P_{3}^{(\text{-}1,\text{-}1)}$:

$\displaystyle
P_{n}^{(\alpha,\beta)}(x)=\frac{(-1)^{n}}{2^{n}n!}\left(1-x\right)^{-\alpha}\left(1+x\right)^{-\beta}\frac{d^{n}}{dx^{n}}\left(\left(1-x\right)^{\alpha+n}\left(1+x\right)^{\beta+n}\right),\hskip
8.53581pt\text{for}\,\,n\geq 2\hskip 0.7113pt.$

As before we take $\alpha,\beta\rightarrow 1^{+}$ in the above to obtain an
explicit formula for $P_{n}^{(\text{-}1,\text{-}1)}$. ∎

In addition to computing these polynomials via a recurrence relation, it is
also possible to represent each $P_{n}^{(\text{-}1,\text{-}1)}$ as the
difference of two (rescaled) Legendre polynomials. Since the Legendre
polynomials already have efficient implementations in the majority of high-
level programming languages, this second approach is particularly appealing.

######  Theorem 12 (Relationship between the Jacobi-like and Legendre
polynomials).

For $n\geq 1$, we have

$\displaystyle P_{n+1}^{(\text{-}1,\text{-}1)}(x)=\frac{n}{4n+2}\hskip
0.7113pt\big{(}Q_{n+1}(x)-Q_{n-1}(x)\big{)},$

where $Q_{k}$ denotes the $k$-th degree Legendre polynomial defined on
$[-1,1]$.

###### Proof.

Recall that $\frac{d}{dx}\big{(}P^{(\text{-}1,\text{-}1)}_{n+1}\hskip
0.7113pt\big{)}=\frac{n}{2}P^{(0,0)}_{n}$ for $n\geq 1$, where
$P^{(0,0)}_{n}(=Q_{n})$ is the $n$-th degree Legendre polynomial. Therefore
differentiating both sides of (20) yields

$\displaystyle\frac{1}{2}\hskip 0.7113ptn(n+1)(n+2)\hskip 0.7113ptQ_{n+1}(x)$
$\displaystyle=(n+1)(2n+1)P_{n+1}^{(\text{-}1,\text{-}1)}(x)+\frac{1}{2}\hskip
0.7113ptn(n+1)(2n+1)\hskip 0.7113ptx\hskip 0.7113ptQ_{n}(x)$
$\displaystyle\hskip 14.22636pt-\frac{1}{2}\hskip 0.7113ptn(n-1)(n+1)\hskip
0.7113ptQ_{n-1}(x)\hskip 0.7113pt.$

Hence by simplifying and rearranging the above, we have that for $n\geq 1$,

$\displaystyle(2n+1)P_{n+1}^{(\text{-}1,\text{-}1)}(x)$
$\displaystyle=\frac{1}{2}\hskip 0.7113ptn\hskip
0.7113pt\big{(}Q_{n+1}(x)-Q_{n-1}(x)\big{)}$ $\displaystyle\hskip
14.22636pt+\frac{1}{2}\hskip 0.7113ptn\hskip 0.7113pt\big{(}(n+1)\hskip
0.7113ptQ_{n+1}(x)-(2n+1)\hskip 0.7113ptx\hskip 0.7113ptQ_{n}(x)+n\hskip
0.7113ptQ_{n-1}(x)\big{)}.$

We see the last term is zero by a recurrence relation for Legendre polynomials
[20]. ∎

In addition to viewing the polynomials $\\{e_{k}\\}$ as orthogonal with
respect to the weight function $w(x):=\frac{1}{x(1-x)}$ , we can characterize
them via their iterated time integrals. In particular, for $1\leq k\leq n-1$ ,
it follows from the integration by parts formula that

$\displaystyle\hskip
28.45274pt\int_{0<s_{1}<\cdots<s_{k}<1}e_{n}(s_{1})\,ds_{1}\cdots\,ds_{k}$
$\displaystyle=\int_{0}^{1}e_{n}(s)\,d\left(\frac{1}{k!}s^{k}\right)$
$\displaystyle=\frac{1}{(k-1)!}\int_{0}^{1}s^{k-1}e_{n}(s)\,ds$
$\displaystyle=-\frac{1}{k!}\int_{0}^{1}s^{k}e_{n}^{\prime}(s)\,ds$
$\displaystyle=0\hskip 0.7113pt,\hskip 55.48286pt(\text{by the orthogonality
of $e_{n}^{\prime}$}).$

Hence for $k\geq 1$, $e_{k}$ is a polynomial with degree $k+1$ that has roots
at $0$ and $1$ as well as $k-1$ trivial iterated integrals against time. By
additionally specifying the $k$-th iterated time integral, it is then possible
to characterize the $k$-th polynomial $e_{k}$.

To conclude this section, we will address the relationship between the
orthogonal Jacobi-like polynomials $\\{e_{k}\\}$ and the Alpert-Rokhlin
wavelets given in definition 4. Since each $e_{k}^{\prime}$ is proportional to
the $k$-th shifted Legendre polynomial, the family of polynomials
$\\{e_{k}^{\prime}\\}$ is orthogonal with respect to the standard
$L^{2}([\hskip 0.28453pt0,1])$ inner product. This orthogonality is exactly
what is needed to satisfy the conditions (4) and (5). Hence for any $q\geq 1$
there exists an Alpert-Rokhlin mother function of order $q$ that is a
piecewise polynomial where both pieces can be rescaled and translated to give
$e_{q-1}^{\prime}$.

##  3 Applications to SDEs

Consider the Stratonovich SDE on the interval $[\hskip 0.28453pt0,T]$

$\displaystyle dy_{t}$ $\displaystyle=f_{0}(y_{t})\,dt+f_{1}(y_{t})\circ
dW_{t}\hskip 0.35565pt,$ (21) $\displaystyle y_{0}$ $\displaystyle=\xi,$

where $\xi\in\mathbb{R}^{e}$ and $f_{i}$ denote bounded $C^{\infty}$ vector
fields on $\mathbb{R}^{e}$ with bounded derivatives. It then follows from the
standard Picard iteration argument that there exists a unique strong solution
$y$ to (21). An important tool in the numerical analysis of this solution is
the stochastic Taylor expansion (see chapter 5 of [22] for a comprehensive
review). For the purposes of this paper, we only require the following
specific Taylor expansion.

######  Theorem 13 (High order Stratonovich-Taylor expansion).

Let $y$ denote the unique strong solution to (21) and let $0\leq s\leq t$.
Then $y_{t}$ can be expanded as follows:

$\displaystyle y_{t}$ $\displaystyle=y_{s}+f_{0}(y_{s})\hskip
0.7113pth+f_{1}(y_{s})\hskip 0.7113ptW_{s,t}+\frac{1}{2}\hskip
0.7113ptf_{1}^{\prime}(y_{s})f_{1}(y_{s})\hskip
0.7113ptW_{s,t}^{2}+\frac{1}{2}\hskip
0.7113ptf_{0}^{\prime}(y_{s})f_{0}(y_{s})\hskip 0.7113pth^{2}$ (22)
$\displaystyle\hskip
9.95845pt+f_{0}^{\prime}(y_{s})f_{1}(y_{s})\int_{s}^{t}\int_{s}^{u}\circ\,dW_{v}\,du+f_{1}^{\prime}(y_{s})f_{0}(y_{s})\int_{s}^{t}\int_{s}^{u}\,dv\circ\,dW_{u}$
$\displaystyle\hskip 9.95845pt+\frac{1}{6}\hskip 0.7113pt\big{(}\hskip
0.7113ptf_{1}^{\prime}(y_{s})f_{1}^{\prime}(y_{s})f_{1}(y_{s})+f_{1}^{\prime\prime}(y_{s})(f_{1}(y_{s}),f_{1}(y_{s}))\big{)}W_{s,t}^{3}$
$\displaystyle\hskip 9.95845pt+\big{(}\hskip
0.7113ptf_{0}^{\prime}(y_{s})f_{1}^{\prime}(y_{s})f_{1}(y_{s})+f_{0}^{\prime\prime}(y_{s})(f_{1}(y_{s}),f_{1}(y_{s}))\big{)}\int_{s}^{t}\int_{s}^{u}\int_{s}^{v}\circ\,dW_{r}\circ\,dW_{v}\,du$
$\displaystyle\hskip 9.95845pt+\big{(}\hskip
0.7113ptf_{1}^{\prime}(y_{s})f_{0}^{\prime}(y_{s})f_{1}(y_{s})+f_{1}^{\prime\prime}(y_{s})(f_{0}(y_{s}),f_{1}(y_{s}))\big{)}\int_{s}^{t}\int_{s}^{u}\int_{s}^{v}\circ\,dW_{r}\,dv\circ
dW_{u}$ $\displaystyle\hskip 9.95845pt+\big{(}\hskip
0.7113ptf_{1}^{\prime}(y_{s})f_{1}^{\prime}(y_{s})f_{0}(y_{s})+f_{1}^{\prime\prime}(y_{s})(f_{1}(y_{s}),f_{0}(y_{s}))\big{)}\int_{s}^{t}\int_{s}^{u}\int_{s}^{v}\,dr\circ
dW_{v}\circ dW_{u}$ $\displaystyle\hskip 9.95845pt+\frac{1}{24}\hskip
0.7113pt\big{(}\hskip
0.7113ptf_{1}^{\prime}(y_{s})f_{1}^{\prime}(y_{s})f_{1}^{\prime}(y_{s})f_{1}(y_{s})+f_{1}^{\prime}(y_{s})f_{1}^{\prime\prime}(y_{s})(f_{1}(y_{s}),f_{1}(y_{s}))$
$\displaystyle\hskip 41.54103pt+\hskip 0.7113pt3\hskip
0.7113ptf_{1}^{\prime\prime}(y_{s})(f_{1}^{\prime}(y_{s})f_{1}(y_{s}),f_{1}(y_{s}))+f_{1}^{\prime\prime\prime}(y_{s})(f_{1}(y_{s}),f_{1}(y_{s}),f_{1}(y_{s}))\big{)}\hskip
0.7113ptW_{s,t}^{4}$ $\displaystyle\hskip 9.95845pt+R(h,y_{s}),$

where $h:=t-s$ and the remainder term has the following uniform estimate for
$h<1$,

$\displaystyle\sup_{y_{s}\in\mathbb{R}^{e}}\left\|R(h,y_{s})\right\|_{L^{2}(\mathbb{P})}\leq
C\,h^{\frac{5}{2}},$ (23)

where the constant $C>0$ depends only on the vector fields of the differential
equation.

From a numerical perspective, the most challenging terms presented in (22) are
those that involve non-trivial third order iterated integrals of Brownian
motion and time. Moreover, the most significant source of discretization error
that high order numerical methods will experience is generally due to
approximating these stochastic integrals. By representing Brownian motion as a
(random) polynomial plus independent noise, we shall derive a new optimal and
unbiased estimator for these third order integrals.

######  Theorem 14.

Let $W$ denote a standard real-valued Brownian motion on $[\hskip
0.28453pt0,1]$. Let $W^{n}$ be the unique $n$-th degree random polynomial with
a root at $0$ and satisfying

$\displaystyle\int_{0}^{1}u^{k}\,dW_{u}^{n}$
$\displaystyle=\int_{0}^{1}u^{k}\,dW_{u}\hskip 0.35565pt,\hskip
8.53581pt\text{for}\hskip 5.69054ptk=0,1,\cdots,n-1.$

Then $W=\hskip 0.7113ptW^{n}+Z^{n}$ , where $Z^{n}$ is a centered Gaussian
process independent of $W^{n}$.

Furthermore, $Z^{n}$ has the following covariance function:

$\displaystyle\operatorname{cov}(Z^{n}_{s},Z^{n}_{t})=K_{B}(s,t)-\sum_{k=1}^{n-1}\lambda_{k}e_{k}(s)e_{k}(t),\hskip
8.53581pt\text{for}\hskip 5.69054pts,t\in[\hskip 0.28453pt0,1],$

where $K_{B}$ denotes the standard Brownian bridge covariance function and
$\\{\lambda_{k}\\}$, $\\{e_{k}\\}$ are the eigenvalues and eigenfunctions that
were defined in the proof of Theorem 5.

###### Proof.

It follows from the integration by parts formula that $W^{n}$ matches the
increment and $n-1$ iterated time integrals of Brownian motion that appear in
(15). Hence $W^{n}$ is also the polynomial defined in Theorem 7 and
$W=W^{n}+Z^{n}$ where

$\displaystyle W^{n}$
$\displaystyle=W_{1}e_{0}+\sum_{k=1}^{n-1}I_{k}e_{k}\hskip 0.35565pt,$
$\displaystyle Z^{n}$ $\displaystyle=\sum_{k=n}^{\infty}I_{k}e_{k}\hskip
0.35565pt.$

Then by Theorem 5, $Z^{n}$ is a centered Gaussian process that is independent
of $W^{n}$. In addition, the covariance function defining $Z^{n}$ can be
directly computed as follows:

$\displaystyle\operatorname{cov}(Z^{n}_{s},Z^{n}_{t})$
$\displaystyle=\operatorname{cov}\left(\sum_{i=n}^{\infty}I_{i}e_{i}(s),\sum_{j=n}^{\infty}I_{j}e_{j}(t)\right)$
$\displaystyle=\sum_{k=n}^{\infty}\lambda_{k}e_{k}(s)e_{k}(t)$
$\displaystyle=K_{B}(s,t)-\sum_{k=1}^{n-1}\lambda_{k}e_{k}(s)e_{k}(t),\hskip
8.53581pt\text{for}\hskip 5.69054pts,t\in[\hskip 0.28453pt0,1].$

Note that the final line is achieved using the representation of $K_{B}$ given
by (8). ∎

The above theorem has an interesting conclusion, namely that there exist
unbiased polynomial approximants of Brownian motion for which the error
process can be independently estimated in an $L^{2}(\mathbb{P})$ sense. In
particular, this theorem already has numerical applications in the case when
$n=2$ and motivates the following definitions:

######  Definition 15.

The standard Brownian parabola $\,\wideparen{W}$ is the unique quadratic
polynomial on $[\hskip 0.28453pt0,1]$ with a root at $0$ and satisfying

$\displaystyle\wideparen{W}_{1}=W_{1}\hskip 0.35565pt,\hskip
14.22636pt\int_{0}^{1}\wideparen{W}_{u}\,du=\int_{0}^{1}W_{u}\,du.$

######  Definition 16.

The standard Brownian arch $\,Z$ is the process $Z:=W-\wideparen{W}$. By
Theorem 14, $Z$ is the centered Gaussian process on $[\hskip 0.28453pt0,1]$
with covariance function

$\displaystyle K_{Z}\left(s,t\right)=\min(s,t)-st-3st(1-s)(1-t),\hskip
8.53581pt\text{for}\hskip 5.69054pts,t\in[\hskip 0.28453pt0,1].$

######  Definition 17.

The rescaled space-time Lévy area of Brownian motion over an interval $[s,t]$
with length $h$ encodes the signed area of the associated bridge process,

$\displaystyle
H_{s,t}:=\frac{1}{h}\int_{s}^{t}W_{s,u}-\frac{u-s}{h}\,W_{s,t}\,du.$

######  Remark 3.18.

Since $e_{1}(t)=\sqrt{6}\,t(1-t)$, we have that $H_{0,1}$ corresponds to
$\frac{\sqrt{6}}{6}I_{1}$ as defined in Theorem 5. Thus,
$H_{s,t}\sim\mathcal{N}\left(0,\frac{1}{12}h\right)$ and $H_{s,t}$ is
independent of $\hskip 0.7113ptW_{s,t}\hskip 0.7113pt$.

By applying the natural scaling of Brownian motion, one can define the
Brownian parabola and Brownian arch processes over any interval $[s,t]$ with
finite size $h=t-s$. Whilst the Brownian arch can be viewed in a similar light
to the Brownian bridge, there are clear qualitative and quantitative
differences in their covariance functions. In particular, the Brownian arch
has less variance at its midpoint compared to most points in $[s,t]$
$\big{(}$by which we mean that
$|\\{u\in[s,t]:\operatorname{Var}(Z_{u})\leq\operatorname{Var}(Z_{\frac{1}{2}(s+t)})\\}|<\frac{1}{2}h\big{)}$.
This is in contrast to the Brownian bridge, which has most variance at its
midpoint. In fact, the Brownian parabola gives a relatively uniform estimate
of the original path.

Fig. 5: Variance profile of the standard Brownian arch.

Using these new definitions, we can study the high order integrals appearing
in (22).

######  Theorem 3.19 (Conditional expectation of a non-trivial Brownian time
integral).

$\displaystyle\mathbb{E}\left[\hskip
0.7113pt\int_{s}^{t}W_{s,u}^{2}\,du\,\Big{|}\,W_{s,t}\hskip
0.35565pt,H_{s,t}\right]=\frac{1}{3}hW_{s,t}^{2}+hW_{s,t}H_{s,t}+\frac{6}{5}hH_{s,t}^{2}+\frac{1}{15}h^{2}.$
(24)

######  Proof 3.20.

By the natural Brownian scaling it is enough to prove the result on $[\hskip
0.28453pt0,1]$. Recall that $W=\wideparen{W}+Z$ where the parabola
$\wideparen{W}$ is completely determined by $\left(W_{1},H_{1}\right)$ and $Z$
is independent of $\left(W_{1},H_{1}\right)$. This leads to a decomposition
for the LHS of (24).

$\displaystyle\mathbb{E}\left[\hskip
0.7113pt\int_{0}^{1}W_{u}^{2}\,du\,\Big{|}\,W_{1},H_{1}\right]$
$\displaystyle=\mathbb{E}\left[\hskip 0.7113pt\int_{0}^{1}\left(\hskip
0.7113pt\wideparen{W}_{u}+Z_{u}\right)^{2}\,du\,\Big{|}\,W_{1},H_{1}\right]$
$\displaystyle=\mathbb{E}\left[\hskip
0.7113pt\int_{0}^{1}\wideparen{W}_{u}^{2}\,du+2\int_{0}^{1}\wideparen{W}_{u}\hskip
0.7113ptZ_{u}\,du+\int_{0}^{1}Z_{u}^{2}\,du\,\Big{|}\,W_{1},H_{1}\right]$
$\displaystyle=\int_{0}^{1}\wideparen{W}_{u}^{2}\,du+2\int_{0}^{1}\wideparen{W}_{u}\,\mathbb{E}\left[Z_{u}\right]\,du+\int_{0}^{1}\mathbb{E}\left[Z_{u}^{2}\hskip
0.7113pt\right]\,du$
$\displaystyle=\int_{0}^{1}\big{(}uW_{1}+6u(1-u)H_{1}\big{)}^{2}du+\int_{0}^{1}u-u^{2}-3u^{2}(1-u)^{2}\,du.$

The result now follows by evaluating the above integrals.

The above theorem has practical applications for SDE simulation as $W_{s,t}$
and $H_{s,t}$ are independent Gaussian random variables and can be easily
generated or approximated. That said, we should first discuss how the iterated
integrals within (22) are connected.

######  Definition 3.21.

The space-space-time Lévy area of Brownian motion over an interval $[s,t]$ is
defined as

$\displaystyle L_{s,t}$
$\displaystyle:=\frac{1}{6}\left(\int_{s}^{t}\int_{s}^{u}\int_{s}^{v}\circ\,dW_{r}\circ\,dW_{v}\,du-2\int_{s}^{t}\int_{s}^{u}\int_{s}^{v}\circ\,dW_{r}\,dv\circ
dW_{u}\right.$ $\displaystyle\hskip
31.29802pt\left.+\int_{s}^{t}\int_{s}^{u}\int_{s}^{v}\,dr\circ dW_{v}\circ
dW_{u}\right).$

We can interpret $L_{s,t}$ as an area between the processes
$\\{W_{s,u}\\}_{u\in[s,t]}$ and $\\{H_{s,u}\\}_{u\in[s,t]}$. Moreover, rough
path theory provides an algebraic structure (called the log-signature) that
relates $(W_{s,t},H_{s,t},L_{s,t})$ to the iterated integrals of space-time
Brownian motion and ultimately to SDE solutions via the log-ODE method (see
[23] for an overview). For our purposes, it is enough to give formulae
relating these Lévy areas to integrals.

######  Theorem 3.22.

Let $H_{s,t}$ and $L_{s,t}$ denote the Lévy areas of Brownian motion given by
definitions 17 and 3.21 respectively. Then the following integral
relationships hold,

$\displaystyle\int_{s}^{t}\int_{s}^{u}\circ\,dW_{v}\,du$
$\displaystyle=\frac{1}{2}hW_{s,t}+hH_{s,t}\hskip 0.35565pt,$
$\displaystyle\int_{s}^{t}\int_{s}^{u}\,dv\circ\,dW_{u}$
$\displaystyle=\frac{1}{2}hW_{s,t}-hH_{s,t}\hskip 0.35565pt,$
$\displaystyle\int_{s}^{t}\int_{s}^{u}\int_{s}^{v}\circ\,dW_{r}\circ\,dW_{v}\,du$
$\displaystyle=\frac{1}{6}hW_{s,t}^{2}+\frac{1}{2}hW_{s,t}H_{s,t}+L_{s,t}\hskip
0.35565pt,$
$\displaystyle\int_{s}^{t}\int_{s}^{u}\int_{s}^{v}\circ\,dW_{r}\,dv\circ
dW_{u}$ $\displaystyle=\frac{1}{6}hW_{s,t}^{2}-2\hskip 0.35565ptL_{s,t}\hskip
0.35565pt,$ $\displaystyle\int_{s}^{t}\int_{s}^{u}\int_{s}^{v}\,dr\circ
dW{v}\circ dW_{u}$
$\displaystyle=\frac{1}{6}hW_{s,t}^{2}-\frac{1}{2}hW_{s,t}H_{s,t}+L_{s,t}\hskip
0.35565pt.$

######  Proof 3.23.

The result follows from numerous applications of integration by parts.

We can now present the new unbiased estimator for third order iterated
integrals of Brownian motion and time. The proposed estimator is fast to
compute and the best $L^{2}(\mathbb{P})$ approximation of these integrals that
is measurable with respect to $\left(W_{s,t},H_{s,t}\right)$.

######  Theorem 3.24 (Conditional moments of Brownian space-space-time Lévy
area).

$\displaystyle\mathbb{E}\left[L_{s,t}\,\big{|}\,W_{s,t}\hskip
0.35565pt,H_{s,t}\right]$ $\displaystyle=\frac{1}{30}\hskip
0.7113pth^{2}+\frac{3}{5}\hskip 0.7113pthH_{s,t}^{2}\,,$ (25)
$\displaystyle\operatorname{Var}\left(L_{s,t}\,\big{|}\,W_{s,t}\hskip
0.35565pt,H_{s,t}\right)$ $\displaystyle=\frac{11}{25200}\hskip
0.7113pth^{4}+h^{3}\left(\frac{1}{720}\hskip
0.7113ptW_{s,t}^{2}+\frac{1}{700}\hskip 0.7113ptH_{s,t}^{2}\right).$ (26)

######  Proof 3.25.

The expectation (25) is simply a consequence of Theorems 3.19 and 3.22.
Without loss of generality, we will consider the above conditional variance on
$[\hskip 0.28453pt0,1]$. Since $\wideparen{W}$ is determined using the
increment $W_{1}$ and space-time Lévy area $H_{1}$, we have

$\displaystyle\operatorname{Var}\left(\hskip
0.7113pt\int_{0}^{1}W_{u}^{2}\,du\,\Big{|}\,W_{1},H_{1}\right)$
$\displaystyle=\operatorname{Var}\left(\hskip
0.7113pt\int_{0}^{1}\wideparen{W}_{u}^{2}\,du+2\int_{0}^{1}\wideparen{W}_{u}\hskip
0.7113ptZ_{u}\,du+\int_{0}^{1}Z_{u}^{2}\,du\,\Big{|}\,W_{1},H_{1}\right)$
$\displaystyle=4\operatorname{Var}\left(\hskip
0.7113pt\int_{0}^{1}\wideparen{W}_{u}\hskip
0.7113ptZ_{u}\,du\,\Big{|}\,W_{1},H_{1}\right)+\operatorname{Var}\left(\hskip
0.7113pt\int_{0}^{1}Z_{u}^{2}\,du\right)$ $\displaystyle\hskip
21.33955pt+4\operatorname{cov}\left(\hskip
0.7113pt\int_{0}^{1}\wideparen{W}_{u}\,Z_{u}\,du\hskip
0.7113pt,\int_{0}^{1}Z_{u}^{2}\,du\,\Big{|}\,W_{1}\hskip
0.35565pt,H_{1}\right).$

Recall $Z=\sum\limits_{k=2}^{\infty}I_{k}\hskip 0.7113pte_{k}$ where
$\left\\{I_{k}\right\\}$ are independent centered Gaussian random variables.

In particular, this means that $Z$ and $-Z$ have the same law. Therefore, we
have that

$\displaystyle\mathbb{E}\left[\hskip
0.7113pt\int_{0}^{1}\wideparen{W}_{u}\hskip
0.7113ptZ_{u}\,du\,\Big{|}\,W_{1}\hskip 0.35565pt,H_{1}\hskip
0.7113pt\right]=0\hskip 0.7113pt,$
$\displaystyle\operatorname{cov}\left(\hskip
0.7113pt\int_{0}^{1}\wideparen{W}_{u}\hskip 0.7113ptZ_{u}\,du\hskip
0.7113pt,\int_{0}^{1}Z_{u}^{2}\,du\,\Big{|}\,W_{1}\hskip
0.35565pt,H_{1}\right)=0\hskip 0.7113pt.$

The remaining two terms were resolved with assistance from Wolfram
Mathematica.

$\displaystyle\operatorname{Var}\left(\hskip
0.7113pt\int_{0}^{1}\wideparen{W}_{u}\hskip
0.7113ptZ_{u}\,du\,\Big{|}\,W_{1}\hskip 0.35565pt,H_{1}\right)$
$\displaystyle\hskip
28.45274pt=\int_{0}^{1}\int_{0}^{1}\wideparen{W}_{u}\hskip
0.7113pt\wideparen{W}_{v}\,\mathbb{E}\left[Z_{u}\hskip
0.7113ptZ_{v}\,|\,W_{1}\hskip 0.35565pt,H_{1}\right]\,du\,dv$
$\displaystyle\hskip
28.45274pt=\int_{0}^{1}\int_{0}^{1}\wideparen{W}_{u}\hskip
0.7113pt\wideparen{W}_{v}\big{(}\min(u,v)-uv-3uv(1-u)(1-v)\big{)}\,du\,dv$
$\displaystyle\hskip
28.45274pt=2\int_{0}^{1}\wideparen{W}_{v}\int_{0}^{v}u\hskip
0.7113pt\wideparen{W}_{u}\,du\,dv-\left(\int_{0}^{1}u\hskip
0.7113pt\wideparen{W}_{u}\,du\right)^{2}-3\left(\hskip
0.7113pt\int_{0}^{1}u(1-u)\wideparen{W}_{u}\,du\right)^{2}$
$\displaystyle\hskip 28.45274pt=2\left(\hskip 0.35565pt\frac{1}{15}\hskip
0.7113ptW_{1}^{2}+\frac{13}{60}\hskip 0.7113ptW_{1}H_{1}+\frac{13}{70}\hskip
0.7113ptH_{1}^{2}\right)-\left(\hskip 0.35565pt\frac{1}{3}\hskip
0.7113ptW_{1}+\frac{1}{2}\hskip 0.7113ptH_{1}\right)^{2}-3\left(\hskip
0.35565pt\frac{1}{12}\hskip 0.7113ptW_{1}+\frac{1}{5}\hskip
0.7113ptH_{1}\right)^{2}$ $\displaystyle\hskip 28.45274pt=\frac{1}{720}\hskip
0.7113ptW_{1}^{2}+\frac{1}{700}\hskip 0.7113ptH_{1}^{2}\hskip 0.7113pt.$
$\displaystyle\operatorname{Var}\left(\int_{0}^{1}Z_{u}^{2}\,du\right)$
$\displaystyle=\mathbb{E}\left[\left(\hskip
0.7113pt\int_{0}^{1}Z_{u}^{2}\,du\right)^{2}\hskip
0.7113pt\right]-\left(\mathbb{E}\left[\hskip
0.7113pt\int_{0}^{1}Z_{u}^{2}\,du\right]\right)^{2}$
$\displaystyle=\int_{0}^{1}\int_{0}^{1}\mathbb{E}\left[Z_{u}^{2}\hskip
0.7113ptZ_{v}^{2}\right]\,du\,dv-\left(\mathbb{E}\left[\hskip
0.7113pt\int_{0}^{1}Z_{u}^{2}\,du\right]\right)^{2}$
$\displaystyle=\int_{0}^{1}\int_{0}^{1}\mathbb{E}\left[Z_{u}^{2}\hskip
0.7113pt\right]\mathbb{E}\left[Z_{v}^{2}\hskip 0.7113pt\right]+2\hskip
0.7113pt\big{(}\hskip 0.7113pt\mathbb{E}[Z_{u}\hskip 0.7113ptZ_{v}]\hskip
0.7113pt\big{)}^{2}\,du\,dv-\left(\mathbb{E}\left[\hskip
0.7113pt\int_{0}^{1}Z_{u}^{2}\,du\right]\right)^{2}$
$\displaystyle=4\int_{0}^{1}\int_{0}^{v}\big{(}u-uv-3uv(1-u)(1-v)\big{)}^{2}\,du\,dv$
$\displaystyle=\frac{11}{6300}\hskip 0.7113pt.$

By Theorem 3.22, the above gives an explicit formula for the conditional
variance (26).

$\displaystyle\operatorname{Var}\left(L_{1}\,\big{|}\,W_{1}\hskip
0.35565pt,H_{1}\right)=\operatorname{Var}\left(\frac{1}{2}\int_{0}^{1}W_{u}^{2}\,du\,\Big{|}\,W_{1}\hskip
0.35565pt,H_{1}\right)=\frac{11}{25200}+\frac{1}{720}\hskip
0.7113ptW_{1}^{2}+\frac{1}{700}\hskip 0.7113ptH_{1}^{2}\hskip 0.7113pt.$

By the natural Brownian scaling, the result on the interval $[s,t]$ directly
follows.

######  Remark 3.26.

The conditional variance (26) allows one to estimate local $L^{2}(\mathbb{P})$
errors for certain numerical methods and thus may be useful when choosing step
sizes.

Therefore in order to propagate a numerical solution of (21) over an interval
$[s,t]$, one can generate $\left(W_{s,t},H_{s,t}\right)$ exactly and then
approximate $L_{s,t}$ using Theorem 3.24. However, there are many numerical
methods that could be used to solve a given SDE.

###  3.1 Examples of ODE methods

We will consider the following two methods:

######  Definition 3.27 (High order log-ODE method).

For a fixed number of steps $N$ we can construct a numerical solution
$\left\\{Y_{k}\right\\}_{0\leq k\leq N}$ of (21) by setting $Y_{0}:=\xi$ and
for each $k\in[0\mathrel{{.}\,{.}}\nobreak N-1]$, defining $Y_{k+1}$ to be the
solution at $u=1$ of the following ODE:

$\displaystyle\frac{dz}{du}$ $\displaystyle=f_{0}(z)\hskip
0.7113pth+f_{1}(z)\hskip 0.7113ptW_{t_{k},t_{k+1}}+[f_{1},f_{0}](z)\cdot
hH_{t_{k},t_{k+1}}$ (27) $\displaystyle\hskip
41.25648pt+[f_{1},[f_{1},f_{0}]](z)\cdot\mathbb{E}\big{[}\hskip
0.7113ptL_{t_{k},t_{k+1}}\,\big{|}\,W_{t_{k},t_{k+1}}\hskip
0.35565pt,H_{t_{k},t_{k+1}}\hskip 0.7113pt\big{]},$ $\displaystyle z_{0}$
$\displaystyle=Y_{k}\hskip 0.35565pt,$

where $h:=\frac{T}{N}$, $t_{k}:=kh$ and $[\,\cdot\hskip 1.42262pt,\cdot\,]$
denotes the standard Lie bracket of vector fields.

######  Definition 3.28 (The parabola-ODE method).

For a fixed number of steps $N$ we can construct a numerical solution
$\left\\{Y_{k}\right\\}_{0\leq k\leq N}$ of (21) by setting $Y_{0}:=\xi$ and
for each $k\in[0\mathrel{{.}\,{.}}\nobreak N-1]$, defining $Y_{k+1}$ to be the
solution at $u=1$ of the following ODE:

$\displaystyle\frac{dz}{du}$ $\displaystyle=f_{0}(z)\hskip
0.7113pth+f_{1}(z)\hskip 0.7113pt\big{(}W_{t_{k},t_{k+1}}+(6-12u)\hskip
0.7113ptH_{t_{k},t_{k+1}}\big{)},$ (28) $\displaystyle z_{0}$
$\displaystyle=Y_{k}\hskip 0.35565pt,$

where $h:=\frac{T}{N}$ and $t_{k}:=kh$.

In both numerical methods the true solution $y$ at time $t_{k}$ can be
approximated by $Y_{k}$. Whilst there are different ways of interpolating
between the successive approximations $Y_{k}$ and $Y_{k+1}$, for this paper we
will simply interpolate between such points linearly. To analyse the above
methods, we shall first note the key differences between them. The first
important distinction between the two methods is a purely practical one.
Although both of these methods involve computing a numerical solution of an
ODE, the parabola method does not require one to explicitly resolve vector
field derivatives. The second significant difference can be seen in the Taylor
expansions of the methods.

######  Theorem 3.29.

Let $Y^{\text{log}}$ be the one-step approximation defined by the log-ODE
method on the interval $[s,t]$ with initial value $Y_{0}^{\text{log}}=y_{s}$.
Then for sufficiently small $h$

$\displaystyle Y_{1}^{\text{log}}=y_{t}-[f_{1},[f_{1},f_{0}]](y_{s})\hskip
1.42262pt\big{(}\hskip 0.7113ptL_{s,t}-\mathbb{E}\big{[}\hskip
0.7113ptL_{s,t}\,\big{|}\,W_{s,t}\hskip 0.35565pt,H_{s,t}\hskip
0.7113pt\big{]}\big{)}+O(h^{\frac{5}{2}}).$ (29)

Similarly, let $Y^{\text{para}}$ denote the one-step approximation given by
the parabola-ODE method on the interval $[s,t]$ with the same initial value.
Then for sufficiently small $h$

$\displaystyle
Y_{1}^{\text{para}}=y_{t}-[f_{1},[f_{1},f_{0}]](y_{s})\left(L_{s,t}-\frac{3}{5}H_{s,t}^{2}\right)+O(h^{\frac{5}{2}}).$
(30)

Note that $O(h^{\frac{5}{2}})$ denotes terms which can be estimated in an
$L^{2}(\mathbb{P})$ sense as in (23).

######  Proof 3.30.

In order to derive (29), we must compute the Taylor expansion of (27). Let $F$
denote the vector field defined in (27) that was constructed from $f_{0}$ and
$f_{1}$. Then $F$ is smooth, and it follows from the classical Taylor’s
theorem for ODEs that

$\displaystyle Y_{1}^{\text{log}}$
$\displaystyle=y_{s}+F(y_{s})+\frac{1}{2}F^{\prime}(y_{s})F(y_{s})+\frac{1}{6}F^{\prime\prime}(y_{s})(F(y_{s}),F(y_{s}))+\frac{1}{6}F^{\prime}(y_{s})F^{\prime}(y_{s})F(y_{s})$
$\displaystyle\hskip
22.76219pt+\frac{1}{24}F^{\prime}(y_{s})F^{\prime}(y_{s})F^{\prime}(y_{s})F(y_{s})+\frac{1}{24}F^{\prime}(y_{s})F^{\prime\prime}(y_{s})(F(y_{s}),F(y_{s}))$
$\displaystyle\hskip
22.76219pt+\,\frac{1}{8}F^{\prime\prime}(y_{s})(F^{\prime}(y_{s})F(y_{s}),F(y_{s}))+\frac{1}{24}F^{\prime\prime\prime}(y_{s})(F(y_{s}),F(y_{s}),F(y_{s}))$
$\displaystyle\hskip
22.76219pt+\frac{1}{24}\int_{0}^{1}(1-u)^{4}\frac{d^{5}}{du^{5}}\big{(}Y^{\text{log}}\big{)}(u)\,du.$

We shall first consider the remainder term, which can be directly estimated as
follows:

$\displaystyle\left\|\,\int_{0}^{1}(1-u)^{4}\frac{d^{5}}{du^{5}}\big{(}Y^{\text{log}}\big{)}(u)\,du\,\right\|_{L^{2}(\mathbb{P})}$
$\displaystyle\leq\sup_{u\in[\hskip
0.28453pt0,1]}\left\|\frac{d^{5}}{du^{5}}\big{(}Y^{\text{log}}(u)\big{)}\right\|_{L^{2}(\mathbb{P})}.$

One can define the degree of each term in the above Taylor expansion by
counting the number of times functions from
$\\{F,F^{\prime},F^{\prime\prime},\cdots\\}$ appear. Therefore, after
expanding the fifth derivative of $Y^{\text{log}}$ we can see that the
remainder term has a degree of five. Since the largest component of $F$ is
$f_{1}(\cdot)\hskip 0.7113ptW_{s,t}$, both $F$ and its derivatives are
$O(h^{\frac{1}{2}})$. Hence the remainder term in the above Taylor expansion
will be $O(h^{\frac{5}{2}})$ as in (23). Moreover, the only terms of degree
four that are not $O(h^{\frac{5}{2}})$ are those involving $W_{s,t}^{4}$.
It is now enough to analyse just the terms appearing in the first line of the
expansion. By substituting the formula for $F$ given by (27) into the first
line and then rearranging the resulting terms, we can obtain a Taylor
expansion for $Y_{1}^{\text{log}}$ that resembles (22) as

$\displaystyle F(y_{s})$ $\displaystyle=f_{0}(y_{s})\hskip
0.7113pth+f_{1}(y_{s})\hskip
0.7113ptW_{s,t}+\big{(}f_{0}^{\prime}(y_{s})f_{1}(y_{s})-f_{1}^{\prime}(y_{s})f_{0}(y_{s})\big{)}\hskip
0.7113pthH_{s,t}$ $\displaystyle\hskip
14.22636pt+\big{(}f_{0}^{\prime}(y_{s})f_{1}^{\prime}(y_{s})f_{1}(y_{s})+f_{0}^{\prime\prime}(y_{s})(f_{1}(y_{s}),f_{1}(y_{s}))$
$\displaystyle\hskip 28.45274pt-2\hskip
0.7113ptf_{1}^{\prime}(y_{s})f_{0}^{\prime}(y_{s})f_{1}(y_{s})-2\hskip
0.7113ptf_{1}^{\prime\prime}(y_{s})(f_{0}(y_{s}),f_{1}(y_{s}))$
$\displaystyle\hskip
28.45274pt+f_{1}^{\prime}(y_{s})f_{1}^{\prime}(y_{s})f_{0}(y_{s})+f_{1}^{\prime\prime}(y_{s})(f_{1}(y_{s}),f_{0}(y_{s}))\big{)}\left(\frac{1}{30}\hskip
0.7113pth^{2}+\frac{3}{5}\hskip 0.7113pthH_{s,t}^{2}\right),$ $\displaystyle
F^{\prime}(y_{s})F(y_{s})$
$\displaystyle=f_{1}^{\prime}(y_{s})f_{1}(y_{s})\hskip
0.7113ptW_{s,t}^{2}+\big{(}f_{0}^{\prime}(y_{s})f_{1}(y_{s})+f_{1}^{\prime}(y_{s})f_{0}(y_{s})\big{)}\hskip
0.7113pthW_{s,t}$ $\displaystyle\hskip
14.22636pt+\big{(}f_{0}^{\prime}(y_{s})f_{1}^{\prime}(y_{s})f_{1}(y_{s})-f_{1}^{\prime}(y_{s})f_{1}^{\prime}(y_{s})f_{0}(y_{s})$
$\displaystyle\hskip
28.45274pt+f_{0}^{\prime\prime}(y_{s})(f_{1}(y_{s}),f_{1}(y_{s}))-f_{1}^{\prime\prime}(y_{s})(f_{0}(y_{s}),f_{1}(y_{s}))\big{)}\hskip
0.7113pthW_{s,t}H_{s,t}$ $\displaystyle\hskip
14.22636pt+f_{0}^{\prime}(y_{s})f_{0}(y_{s})\hskip
0.7113pth^{2}+O\big{(}h^{\frac{5}{2}}\big{)},$ $\displaystyle
F^{\prime}(y_{s})F^{\prime}(y_{s})F(y_{s})$
$\displaystyle=f_{1}^{\prime}(y_{s})f_{1}^{\prime}(y_{s})f_{1}(y_{s})\hskip
0.7113ptW_{s,t}^{3}+f_{0}^{\prime}(y_{s})f_{1}^{\prime}(y_{s})f_{1}(y_{s})\hskip
0.7113pthW_{s,t}^{2}$ $\displaystyle\hskip
14.22636pt+\big{(}f_{1}^{\prime}(y_{s})f_{0}^{\prime}(y_{s})f_{1}(y_{s})+f_{1}^{\prime}(y_{s})f_{1}^{\prime}(y_{s})f_{0}(y_{s})\big{)}\hskip
0.7113pthW_{s,t}^{2}+O\big{(}h^{\frac{5}{2}}\big{)},$ $\displaystyle
F^{\prime\prime}(y_{s})\big{(}F(y_{s}),F(y_{s})\big{)}$
$\displaystyle=f_{1}^{\prime\prime}(y_{s})(f_{1}(y_{s}),f_{1}(y_{s}))\hskip
0.7113ptW_{s,t}^{3}+2\hskip
0.7113ptf_{1}^{\prime\prime}(y_{s})(f_{0}(y_{s}),f_{1}(y_{s}))\hskip
0.7113pthW_{s,t}^{2}$ $\displaystyle\hskip
14.22636pt+f_{0}^{\prime\prime}(y_{s})(f_{1}(y_{s}),f_{1}(y_{s}))\hskip
0.7113pthW_{s,t}^{2}+O\big{(}h^{\frac{5}{2}}\big{)}.$

Therefore, by summing the above formulae for $F$ (and its derivatives) we can
derive an expansion of $Y_{1}^{\text{log}}$ in terms of $f_{0},f_{1}$ and
$\big{(}h,W_{s,t},H_{s,t}\big{)}$ that has an $O\big{(}h^{\frac{5}{2}}\big{)}$
remainder. By comparing this with the stochastic Taylor expansion (22), the
result (29) follows.

Arguing (30) is fairly straightforward and does not require extensive
computations. Using the substitution $\wideparen{Y}_{u}=z_{\frac{1}{h}(u-s)}$
for $u\in[s,t]$, the ODE (28) can be rewritten as

$\displaystyle d\wideparen{Y}_{u}$ $\displaystyle=f_{0}\big{(}\hskip
1.42262pt\wideparen{Y}_{u}\big{)}\hskip 0.7113ptdu+f_{1}\big{(}\hskip
1.42262pt\wideparen{Y}_{u}\big{)}\hskip 0.7113ptd\wideparen{W}_{u}\hskip
0.35565pt,$ (31) $\displaystyle\wideparen{Y}_{s}$ $\displaystyle=y_{s}\hskip
0.35565pt,$

where $\wideparen{W}$ denotes the Brownian parabola defined by
$(W_{s,t},H_{s,t})$ on the interval $[s,t]$.

By emulating the derivation of the Stratonovich-Taylor expansion (22), it is
possible to Taylor expand (31) in the same fashion. The only difference is
that Stratonovich integrals with respect to $W$ are replaced with Riemann-
Stieltjes integrals against $\wideparen{W}$.

In particular, by the change-of-variable formula for ODEs (exercise 3.17)
given in [24], we see that the remainder term of such a Taylor expansion will
have the below form:

$\displaystyle\wideparen{R}\,=\sum_{\begin{subarray}{c}i_{1},\hskip
0.7113pt\cdots,\hskip 0.7113pti_{n}\in\\{0,1\\}\\\ i_{1}\hskip 0.7113pt+\hskip
0.7113pt\cdots\hskip 0.7113pt+\hskip 0.7113pti_{n}\hskip 0.7113pt=\hskip
0.99585pt2n-4\end{subarray}}\int_{s\hskip 0.7113pt<\hskip 0.35565ptr_{1}\hskip
0.35565pt<\hskip 0.35565pt\cdots\hskip 0.35565pt<\hskip 0.35565ptr_{n}<\hskip
0.7113ptt}f_{i_{1},\hskip 0.7113pt\cdots,\hskip 0.7113pti_{n}}\big{(}\hskip
1.42262pt\wideparen{Y}_{r_{1}}\big{)}-f_{i_{1},\hskip 0.7113pt\cdots,\hskip
0.7113pti_{n}}\big{(}\hskip 1.42262pt\wideparen{Y}_{s}\hskip
0.7113pt\big{)}\,d\wideparen{W}_{r_{1}}^{i_{1}}\,\cdots\,d\wideparen{W}_{r_{n}}^{i_{n}}\hskip
0.7113pt,$

where we have identified an additional “ zero” coordinate of $\wideparen{W}$
with time, $\wideparen{W}_{t}^{0}:=t$, and for each index $(i_{1},\cdots
i_{n})$, the function $f_{i_{1},\hskip 0.7113pt\cdots,\hskip
0.7113pti_{n}}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ consists of finitely
many compositions of $f_{0},f_{1}$ along with their derivatives (and thus is
Lipschitz continuous).

Therefore each term in the expansion of (31) can be estimated in
$L^{2}(\mathbb{P})$ by applying the natural Brownian scaling to the
corresponding iterated integral of $\wideparen{W}$ with time. As before, the
largest differences are the $O(h^{2})$ terms involving third order integrals.
Fortunately, iterated integrals of the Brownian parabola can be computed
explicitly:

$\displaystyle\int_{s}^{t}\int_{s}^{u}\int_{s}^{v}\circ\,dW_{r}\circ\,dW_{v}\,du-\int_{s}^{t}\int_{s}^{u}\int_{s}^{v}\,d\wideparen{W}_{r}\,d\wideparen{W}_{v}\,du$
$\displaystyle=L_{s,t}-\frac{3}{5}H_{s,t}^{2}\hskip 0.35565pt.$

The result (30) is now a direct consequence of Theorem 3.22 along with the
above.

Theorem 3.29 shows that both methods give a one-step approximation error of
$O(h^{2})$. This means that the log-ODE and parabola-ODE methods are both
locally high order; however there is a significant difference in how these
methods propagate local errors. The reason is that the $O(h^{2})$ components
of the log-ODE local errors give a martingale, whilst the $O(h^{2})$ part for
each parabola-ODE local error has non-zero expectation. Thus the log-ODE
method is globally high order whilst the parabola method is not. However,
since the parabola-ODE method is straightforward to implement and locally high
order, one could expect it to perform well compared to other low order
methods. In the numerical example, we shall see that the parabola method has
the same order of convergence as the piecewise linear approach but gives
significantly smaller errors. To conclude this section, we will present the
orders of convergence for both methods.

######  Definition 3.31 (Strong convergence).

A numerical solution $Y$ for (21) is said to converge in a strong sense with
order $\alpha$ if there exists a constant $C>0$ such that

$\displaystyle\left\|Y_{N}-y_{T}\right\|_{L^{2}\left(\mathbb{P}\right)}\leq
Ch^{\alpha},$

for all sufficiently small step sizes $h=\frac{T}{N}$.

######  Definition 3.32 (Weak convergence).

A numerical solution $Y$ for (21) is said to converge in a weak sense with
order $\beta$ if for any polynomial $p$ there exists $C_{p}>0$ such that

$\displaystyle\big{|}\,\mathbb{E}\left[\,p\hskip 0.7113pt(\hskip
0.7113ptY_{N})\right]-\mathbb{E}\left[\,p\hskip 0.7113pt(\hskip
0.7113pty_{T})\right]\big{|}\leq C_{p}\,h^{\beta},$

for all sufficiently small step sizes $h=\frac{T}{N}$.

######  Theorem 3.33 (Orders of convergence).

For a general SDE (21), the log-ODE method converges in a strong sense with
order 1.5 and a weak sense with order 2.0. The parabola-ODE method converges
in both a strong and weak sense with order 1.0.

######  Proof 3.34.

Note that Theorem 3.29 establishes the Taylor expansions of both methods. The
strong convergence can then be shown as in the proof of Theorem 11.5.1 in
[22]. Moreover, the proof of Theorem 11.5.1 also provides the orders of strong
convergence. Similarly weak convergence follows directly from the Taylor
expansions (29) & (30), and the rate of convergence can be shown as in the
proof of Theorem 14.5.2 in [22].

##  4 A numerical example

We shall demonstrate the ideas presented so far using various discretizations
of Inhomogeneous Geometric Brownian Motion (IGBM)

$\displaystyle dy_{t}=a(b-y_{t})\,dt+\sigma\hskip 0.35565pty_{t}\,dW_{t}\hskip
0.35565pt,$ (32)

where $a\geq 0$ and $b\in\mathbb{R}$ are the mean reversion parameters and
$\sigma\geq 0$ is the volatility. As the vector fields are smooth, the SDE
(32) can be expressed in Stratonovich form:

$\displaystyle dy_{t}=\tilde{a}(\tilde{b}-y_{t})\,dt+\sigma\hskip
0.35565pty_{t}\circ dW_{t}\hskip 0.35565pt,$ (33)

where $\tilde{a}:=a+\frac{1}{2}\sigma^{2}$ and
$\tilde{b}:=\frac{2ab}{2a+\sigma^{2}}$ denote the “adjusted” mean reversion
parameters.

IGBM is an example of a one-factor short rate model and has seen recent
attention in the mathematical finance literature as an alternative to popular
models [16, 17]. IGBM is also one of the simplest SDEs that has no known
method of exact simulation. We will investigate the strong and weak
convergence rates of the following methods:

  1. 1.

Log-ODE method (see definition 3.27)
Since the vector fields of (33) give constant Lie brackets, this method
becomes

$\displaystyle Y_{k+1}^{\text{log}}$
$\displaystyle:=Y_{k}^{\text{log}}e^{-\tilde{a}h+\sigma W_{t_{k},t_{k+1}}}$
$\displaystyle\hskip 11.38109pt+abh\left(1-\sigma
H_{t_{k},t_{k+1}}+\sigma^{2}\left(\frac{3}{5}H_{t_{k},t_{k+1}}^{2}+\frac{1}{30}h\right)\right)\frac{e^{-\tilde{a}h+\sigma
W_{t_{k},t_{k+1}}}-1}{-\tilde{a}h+\sigma W_{t_{k},t_{k+1}}}\,,$ $\displaystyle
Y_{0}^{\text{log}}$ $\displaystyle:=y_{0}\hskip 0.35565pt.$

  2. 2.

Parabola-ODE method (see definition 3.28)
As the SDE (33) is quite analytically tractable, this method is expressible as

$\displaystyle Y_{k+1}^{\text{para}}$ $\displaystyle:=e^{-\tilde{a}h+\sigma
W_{t_{k},t_{k+1}}}\left(Y_{k}^{\text{para}}+ab\int_{t_{k}}^{t_{k+1}}e^{\tilde{a}\left(s-t_{k}\right)-\sigma\wideparen{W}_{t_{k},s}}\,ds\right),$
$\displaystyle Y_{0}^{\text{para}}$ $\displaystyle:=y_{0}\hskip 0.35565pt.$

The integral above will be computed by $3$-point Gauss-Legendre quadrature.

  3. 3.

Piecewise linear method (see [25] for definition and proof of convergence)
Just as above, this method can be simplified to give a straightforward
formula.

$\displaystyle Y_{k+1}^{\text{lin}}$
$\displaystyle:=Y_{k}^{\text{lin}}e^{-\tilde{a}h+\sigma
W_{t_{k},t_{k+1}}}+abh\,\frac{e^{-\tilde{a}h+\sigma
W_{t_{k},t_{k+1}}}-1}{-\tilde{a}h+\sigma W_{t_{k},t_{k+1}}}\,,$ $\displaystyle
Y_{0}^{\text{lin}}$ $\displaystyle:=y_{0}\hskip 0.35565pt.$

  4. 4.

Milstein method (see section 6 of [2] and section 10.3 of [22] for overviews)
For this method, we shall take the positive part to guarantee non-negativity.

$\displaystyle Y_{k+1}^{\text{mil}}$
$\displaystyle:=\Big{(}Y_{k}^{\text{mil}}+\tilde{a}(\tilde{b}-Y_{k}^{\text{mil}})h+\sigma\,Y_{k}^{\text{mil}}\,W_{t_{k},t_{k+1}}+\frac{1}{2}\sigma^{2}\,Y_{k}^{\text{mil}}\,W_{t_{k},t_{k+1}}^{2}\Big{)}^{+},$
$\displaystyle Y_{0}^{\text{mil}}$ $\displaystyle:=y_{0}\hskip 0.35565pt.$

  5. 5.

Euler-Maruyama method (see sections 4, 5 of [2] and section 10.2 of [22])
Just as above, we take the positive part of each step to ensure non-
negativity.

$\displaystyle Y_{k+1}^{\text{eul}}$
$\displaystyle:=\Big{(}Y_{k}^{\text{eul}}+a(b-Y_{k}^{\text{eul}})h+\sigma\,Y_{k}^{\text{eul}}\,W_{t_{k},t_{k+1}}\Big{)}^{+},$
$\displaystyle Y_{0}^{\text{eul}}$ $\displaystyle:=y_{0}\hskip 0.35565pt.$

Note that the explicit formula for the log-ODE method comes from the Lie
brackets:

$\displaystyle[f_{1},f_{0}](y)$
$\displaystyle=f_{0}^{\prime}(y)f_{1}(y)-f_{1}^{\prime}(y)f_{0}(y)$
$\displaystyle=-\tilde{a}\sigma y-\tilde{a}\sigma(\tilde{b}-y)$
$\displaystyle=-ab\sigma\hskip 0.7113pt,$
$\displaystyle[f_{1},[f_{1},f_{0}]](y)$
$\displaystyle=f_{0}^{\prime}(y)f_{1}^{\prime}(y)f_{1}(y)-2f_{1}^{\prime}(y)f_{0}^{\prime}(y)f_{1}(y)+f_{1}^{\prime}(y)f_{1}^{\prime}(y)f_{0}(y)$
$\displaystyle\hskip
22.76219pt+f_{0}^{\prime\prime}(y)(f_{1}(y),f_{1}(y))-2\hskip
0.7113ptf_{1}^{\prime\prime}(y)(f_{0}(y),f_{1}(y))+f_{1}^{\prime\prime}(y)(f_{1}(y),f_{0}(y))$
$\displaystyle=-\tilde{a}\sigma^{2}y+2\hskip
0.7113pt\tilde{a}\sigma^{2}y+\tilde{a}\sigma^{2}(\tilde{b}-y)$
$\displaystyle=ab\sigma^{2}\hskip 0.7113pt,$

and the formula for the parabola-ODE method was derived using a change of
variable.

The Euler-Maruyama and Milstein methods are included in the numerical
experiment as benchmarks to test how the proposed methods compare to well-
known methods. As before, we will be discretizing the SDE over a uniform
partition with mesh size $h$.

Fig. 6: Log-ODE sample paths of IGBM where $a=0.1$, $b=0.04$, $\sigma=0.6$ and
$y_{0}=0.06$.
The above sample paths experience larger fluctuations the further away from
zero they are.

Below is the definition of the error estimators used to analyse the numerical
methods.

######  Definition 4.35 (Strong and weak error estimators).

For each $N\geq 1$, let $Y_{N}$ denote a numerical solution of (32) computed
at time $T$ using a fixed step size $h=\frac{T}{N}$. We can define the
following estimators for quantifying strong and weak convergence:

$\displaystyle S_{N}$ $\displaystyle:=\sqrt{\hskip
0.7113pt\mathbb{E}\left[\left(Y_{N}-Y_{T}^{fine}\right)^{2}\hskip
1.42262pt\right]}\hskip 0.7113pt,$ (34) $\displaystyle E_{N}$
$\displaystyle:=\left|\,\mathbb{E}\big{[}\big{(}Y_{N}-b\hskip
0.7113pt\big{)}^{+}\hskip
0.7113pt\big{]}-\mathbb{E}\big{[}\big{(}Y_{T}^{fine}-b\hskip
0.7113pt\big{)}^{+}\hskip 0.7113pt\big{]}\right|,$ (35)

where the above expectations are approximated by standard Monte-Carlo
simulation and $Y_{T}^{fine}$ is the numerical solution of (32) obtained at
time $T$ using the log-ODE method with a “ fine” step size of
$\min\left(\frac{h}{10},\frac{T}{1000}\right)$. The fine step size is chosen
so that the $L^{2}(\mathbb{P})$ error between $Y_{T}^{fine}$ and the true
solution $y$ is negligible compared to $S_{N}$. Note that $Y_{N}$ and
$Y_{T}^{fine}$ are both computed with respect to the same Brownian paths.

In this numerical example, we shall use the same parameter values as in [16],
namely $a=0.1$, $b=0.04$, $\sigma=0.6$ and $y_{0}=0.06$. We will also fix the
time horizon at $T=5$.

We will now present our results for the numerical experiment that is described
above. (Code for this example can be found at github.com/james-m-foster/igbm-
simulation)

Fig. 7: $S_{N}$ computed with 100,000 sample paths as a function of step size
$h=\frac{T}{N}$.

From the above graph we see that the log-ODE method is by far the most
accurate. This is epitomized by the fact that the numerical error produced by
$100$ steps of the log-ODE method is comparable to the error of the parabola
method with $1000$ steps. In addition, whilst there are three methods that
share the same order of convergence it is evident there are magnitudes of
difference between their respective accuracies. For example, the parabola
method is seven times more accurate than piecewise linear. As one might
expect, the Euler-Maruyama and Milstein schemes both perform poorly.

Nevertheless, in order to truly measure the performance of these numerical
methods, we should consider the computational costs required for achieving a
specified accuracy.

Table 2:  Estimated simulation times for computing 100,000 sample paths that
achieve a
given accuracy using a single-threaded C++ program on a desktop computer.

| Log-ODE | Parabola |  Linear | Milstein | Euler
---|---|---|---|---|---
Estimated time to achieve | 0.179 | 0.405 | 1.47 | 15.4 | 0.437
an accuracy of $S_{N}=10^{-4}$ | (s) | (s) | (s) | (s) | (days)
Estimated time to achieve | 0.827 | 3.90 | 14.9 | 157 | 61.2
an accuracy of $S_{N}=10^{-5}$ | (s) | (s) | (s) | (s) | (days)

The above times are estimated by the graph shown in Fig 7 and the following
table:

Table 3:  Simulation times for computing 100,000 sample paths with $100$ steps
per path
using a single-threaded C++ program on a desktop computer.

| Log-ODE | Parabola |  Linear | Milstein |  Euler
---|---|---|---|---|---
Computation time (s) | 2.44 | 2.95 | 1.48 | 1.18 | 1.17

Finally, we will investigate the rates of weak convergence for these numerical
methods.

Fig. 8: $E_{N}$ computed with 500,000 sample paths as a function of step size
$h=\frac{T}{N}$.

The above graph demonstrates that the log-ODE method is especially well-suited
for weak approximation as it achieves a second order convergence rate in this
example. Surprisingly, the middle three methods exhibit almost identical
convergence profiles. As before, we can estimate the computational time needed
to achieve given accuracies.

Table 4:  Estimated simulation times for computing 100,000 sample paths that
achieve a
given accuracy using a single-threaded C++ program on a desktop computer.
| Log-ODE | Parabola |  Linear | Milstein | Euler
---|---|---|---|---|---
Estimated time to achieve |  $<$ 0.240 | 1.69 | 2.15 | 2.78 | 25.5
an accuracy of $E_{N}=10^{-5}$ | (s) | (s) | (s) | (s) | (s)
Estimated time to achieve | 0.240 | 16.9 | 21.6 | 24.1 | 252
an accuracy of $E_{N}=10^{-6}$ | (s) | (s) | (s) | (s) | (s)

We expect the log-ODE and parabola methods to have about twice the
computational cost as the other methods because each step requires generating
two random variables. Table 3 confirms this and thus sampling may be a
bottleneck for these methods. So overall, the numerical evidence supports our
claim that the high order log-ODE method is currently a state-of-the-art
method for the pathwise discretization of IGBM.

##  5 Conclusion

There are primarily three new results established in this paper:

  * $\bullet$

An efficient strong polynomial approximation of Brownian motion

The main result allows one to construct a “smoother” Brownian motion as a
finite sum of $(\text{-}1,\text{-}1)$-Jacobi polynomials with independent
Gaussian weights. Moreover, it was shown that the approximation is optimal in
a weighted $L^{2}(\mathbb{P})$ sense and the surrounding noise is an
independent centered Gaussian process.

  * $\bullet$

Unbiased approximation of third order Brownian iterated integrals

Iterated integrals of Brownian motion and time are important objects in the
study of SDEs as they appear naturally within stochastic Taylor expansions. We
have derived the $L^{2}(\mathbb{P})$-optimal estimator for a class of such
integrals that is measurable with respect to the path’s increment and space-
time Lévy area.

  * $\bullet$

Simulation of Inhomogeneous Geometric Brownian Motion (IGBM)

IGBM is a mean-reverting short rate model used in mathematical finance and
also one of the simplest SDEs that has no known method of exact simulation.
By incorporating the new iterated integral estimator into the log-ODE method
we have developed a high order state-of-the-art numerical method for IGBM.

Furthermore, the results of this paper naturally lead to the following open
questions:

  * $\bullet$

Which weight functions give “ explicit eigenfunctions” for Brownian motion?
(For example, we could try $w(x)=x$ or $w(x)=\frac{1}{x}$ with
$K_{W}(s,t)=\min(s,t)$)

  * $\bullet$

Is it possible to generalize the main theorem to fractional Brownian motion?

  * $\bullet$

What are the most efficient Runge-Kutta methods for general one-dimensional
SDEs that correctly use the new estimator for third order iterated integrals?

  * $\bullet$

Is this polynomial expansion optimal for approximating Lévy area? (see [13])

  * $\bullet$

Which conditional moments can be computed for a given stochastic integral?

  * $\bullet$

How might we construct a piecewise linear path $\gamma$ with the below
properties?

$\displaystyle 1.\,\,\,\,\gamma_{s}=W_{s}\hskip 0.7113pt,\hskip
4.97922pt\gamma_{t}=W_{t}\hskip 0.7113pt.$ $\displaystyle
2.\,\,\,\int_{s}^{t}\gamma_{s,u}\,du=\int_{s}^{t}W_{s,u}\,du\hskip 0.7113pt.$
$\displaystyle 3.\,\,\,\int_{s}^{t}\gamma_{s,u}^{2}\,du=\mathbb{E}\left[\hskip
0.7113pt\int_{s}^{t}W_{s,u}^{2}\,du\,\Big{|}\,\cdots\right].$

  * $\bullet$

Would this method of construction lead to effective cubature paths? (see [26])

Given such a path, we can approximate (21) with a “ piecewise linear” ODE.

$\displaystyle\frac{dY}{du}=f_{0}(\hskip 0.35565ptY)+f_{1}(\hskip
0.35565ptY)\hskip 1.42262pt\frac{d\gamma}{du}\hskip 1.42262pt.$ (36)

(Along each piece of $\gamma$, we would discretize (36) using an appropriate
solver)

  * $\bullet$

How effective is the above piecewise linear ODE method for simulating SDEs?

  * $\bullet$

Can we extend the approximations given in this paper to the SPDE setting?

## References

  * [1] R. M. Mazo, Brownian Motion: Fluctuations, Dynamics and Applications, Clarendron Press, Oxford, 2002.
  * [2] D. J. Higham, An algorithmic introduction to numerical simulation of stochastic differential equations, SIAM Review, Volume 43, 2001.
  * [3] J. M. C. Clark and R. J. Cameron, The maximum rate of convergence of discrete approximations for stochastic differential equations, Stochastic Differential Systems Filtering and Control, Volume 25 in Lecture Notes in Control and Information Sciences, Springer, 1980.
  * [4] D. S. Grebenkov, D. Belyaev and P. W. Jones, A multiscale guide to Brownian motion, Journal of Physics A, Volume 49, 2015.
  * [5] K. Habermann, A semicircle law and decorrelation phenomena for iterated Kolmogorov loops, https://arxiv.org/abs/1904.11484, 2019.
  * [6] F. Castell and J. G. Gaines, The ordinary differential equation approach to asymptotically efficient schemes for solution of stochastic differential equations, Annales de l’Institut Henri Poincaré, Volume 32, 1996.
  * [7] S. J. A. Malham and A. Wiese, Stochastic Lie group integrators, Siam Journal on Scientific Computing, Volume 30, 2008.
  * [8] J. M. C. Clark, An efficient approximation for a class of stochastic differential equations, Advances in Filtering and Optimal Stochastic Control, Proceedings of the IFIP-WG7/1 Working Conference. Cocoyoc, Mexico, 1982. Lecture Notes in Control and Information Sciences, Volume 42, Berlin, 1982.
  * [9] N. J. Newton, An asymptotically efficient difference formula for solving stochastic differential equations, Stochastics, Volume 19, 1986.
  * [10] N. J. Newton, Asymptotically efficient Runge-Kutta methods for a class of Itô and Stratonovich equations, SIAM Journal on Applied Mathematics, Volume 51, 1991.
  * [11] X. Tang and A. Xiao, Asymptotically optimal approximation of some stochastic integrals and its applications to the strong second-order methods, Advances in Computational Mathematics, Volume 45, 2019.
  * [12] P. Przybylowicz, Linear information for approximation of the Itô integrals, Numerical Algorithms, Volume 52, 2009.
  * [13] A. S. Dickinson, Optimal Approximation of the Second Iterated Integral of Brownian Motion, Stochastic Analysis and Applications, Volume 25, 2007.
  * [14] D. Lee and G. W. Wasilkowski, Approximation of linear functionals on a Banach space with a Gaussian measure, Journal of Complexity, Volume 2, 1986.
  * [15] J. Creutzig, T. Müller-Gronbach and K. Ritter, Free-knot spline approximation of stochastic processes, Journal of Complexity, Volume 23, 2007.
  * [16] L. Capriotti, Y. Jiang and G. Shaimerdenova, Approximation Methods for Inhomogeneous Geometric Brownian Motion, International Journal of Theoretical and Applied Finance, Volume 22, 2019.
  * [17] B. Zhao, C. Yan and S. Hodges, Three One-Factor Processes for Option Pricing with a Mean-Reverting Underlying: The Case of VIX, Financial Review, Volume 54, 2019.
  * [18] G. Beylkin, R. Coifman and V. Rokhlin, Wavelets in Numerical Analysis, Wavelets and Their Applications, Jones and Bartlett, 1992.
  * [19] C. Carmeli, A. De Vito and E. Toigo, Vector valued reproducing kernel Hilbert spaces of integrable functions and Mercer theorem, Analysis and Applications, Volume 4, 2006.
  * [20] B. G. S. Doman, The Classical Orthogonal Polynomials, World Scientific, 2016.
  * [21] E. W. Weisstein, Jacobi Polynomial, From MathWorld - A Wolfram Web Resource, https://mathworld.wolfram.com/JacobiPolynomial.html.
  * [22] P. E. Kloeden and E. Platen, Numerical Solution of Stochastic Differential Equations, Springer, 1992.
  * [23] T. Lyons, Rough paths, signatures and the modelling of functions on streams, Proceedings of the International Congress of Mathematicians, Volume 4, 2014.
  * [24] P. K. Friz and N. B. Victoir, Multidimensional stochastic processes as rough paths: theory and applications, Cambridge University Press, 2010.
  * [25] E. Wong and M. Zakai, On the Convergence of Ordinary Integrals to Stochastic Integrals, Annals of Mathematical Statistics, Volume 36, 1965.
  * [26] T. Lyons and N. Victoir, Cubature on Wiener space, Proceedings of the Royal Society of London Series A: Mathematical, Physical and Engineering Sciences, Volume 460, 2004.