****Keywords**** Free Stochastic Differential Equations, Free Probability Theory, Milstein method, Random Matrix Theory, Stochastic Differential Equations, Strong convergence

****AMS Codes**** 46L53, 46L54, 60H10, 65C30

1 Introduction↩︎

Free stochastic differential equations (fSDE) emerged up based on the free probability theory, developed by D. Voiculescu in the time period of early 1980s to 1990s. This allowed several researchers to show, that in principle the Doeblin-Itô-calculus and tools from random matrix theory (RMT) can be transferred to these non-commutative setting in an appropriate way, see e.g. [1], [2], [3], [4], [5], [6] and references therein. Nevertheless there are certain differences which indicate that a merely literal translation of the classical stochastic calculus hits its limits. While the classical stochastic differential equations are driven by at most vector-valued stochastic processes (e.g. Brownian motion or Lévy processes), here, the state space is an abstract von Neumann algebra with unital, normal, faithful trace. The driving process is a so-called free Brownian motion with values in an abstract finite von Neumann algebra. Thus, we have to encounter the non-commutativity. To get an idea, one should consider the von Neumann algebra of \(N\times N\)-matrices \(M_N(\mathbb{C})\) equipped with the normalized trace \(\frac{1}{N}\mathbb{E}(\text{tr}(\cdot))\). Let \(\left(W_1^{(N)},\dots, W_d^{(N)}\right)\) be an \(d\)-tuple of independent, standard \(M_N(\mathbb{C})\)-valued Brownian motions. Then there is a von Neumann algebra \(\mathcal{A}\), a trace \(\varphi:\mathcal{A}\rightarrow \mathbb{C}\) and \(W_1,\dots,W_d\) freely independent, \(\mathcal{A}\)-valued, self-adjoint processes on \([0,\infty)\) (called semicircular processes) such that \[\varphi\left(P\left(W_{i_1}^{(N)}(t_1),\dots,W_{i_r}^{(N)}(t_r)\right)\right)\rightarrow \varphi\left(P(\left(W_{i_1}(t_1),\dots,W_{i_r}(t_r)\right))\right)\] almost sure as \(N\rightarrow\infty\) for all \(i_1,\dots,i_r\in\{1,\dots,d\}\), \(t_1,\dots t_r \geq 0\) and polynomials \(P\) in \(r\) non-commuting variables. An introduction into free probability and the relation to random matrices can be found e.g. in [7], [8], [2]. P. Biane showed in [5], that a certain scaling of the \(N\times N\) Hermitian matrix Brownian motion converges as \(N\rightarrow \infty\) to a so-called free Brownian motion. This motivates the viewpoint, that often formulas of \(W_1,\dots,W_d\) in an abstract setting can be studied by considering \(W_1^{(N)},\dots, W_n^{(N)}\) and then taking \(N\rightarrow \infty\). It was then P. Biane and R. Speicher who developed ([4], [3], [5]) the framework of stochastic calculus with respect to free Brownian motion. In their introduction they motivated stochastic free calculus from the matrix level with successively taking limits \(N\rightarrow\infty\) (although the paper was developed by a different approach). In [9] the author gives a rigorous matrix-valued stochastic calculus and an Itô formula for \(C^2\) scalar functions of Hermitian matrix-valued Itô-processes.
To the best knowledge of the authors, free stochastic differential equations appeared first in [1] driven by the idea, to obtain new processes from given ones. P. Biane and R. Speicher studied in [5] diffusion equations driven by Brownian motion on matrix level and in the limit \(N\rightarrow\infty\). The free analog of the central limit theorem (see e.g. [8], [10]) then shows, that the free stochastic differential equations are a good approximation and helpful modeling tool for the wide-spread used random matrices. A general study on fSDEs was gained by V. Kargin (kargin?). He approached fSDEs by the Cauchy transform of the solution, obtaining partial differential equations in terms of the resolvent of the solution of the fSDE. Additionally, a Picard-Lindelöf-type existence result was given.
We consider free stochastic differential equations equipped with a single free Brownian motion \((W_t)_{t\geq 0}\), i.e. \[\label{eq:fSDEinIntro} dX_t=a(X_t)dt + \sum_{i=1}^db^i(X_t)dW_tc^i(X_t),\tag{1}\] where \(X_t\in\mathcal{A}\) is a self-adjoint operator and \(a,b^i, c^i\) are operator valued, adapted and continuous (in operator norm) functions of \(X_t\). Examples are the free analog of the Ornstein-Uhlenbeck process \(dX_t=\lambda X_t dt + \sigma dW_t\), the so-called geometric Brownian motion I, \(dX_t=\theta X_t dt + \sqrt{X_t}dW_t\sqrt{X_t}\), the Geometric Brownian motion II \(dX_t=\theta X_t dt + X_tdW_t + dW_tX_t\) (we refer for details to kargin?). Note that formulating the process Geometric Brownian motion II via 1 we set \(d=2\), \(b^1(X_t)=c^2(X_t)=X_t\) and \(b^2(X_t)=c^1(X_t)=\text{id}\). For details about the formalism, we refer to 2.3.
The Ornstein-Uhlenbeck free process was also studied in [11]. Free Wishart processes were discussed in [12]. A fSDE appeared in the context of free Jacobi processes in [13]. In [14], the authors stated a free variant of Cox–Ingersoll–Ross (CIR) process. In [15] fSDEs appeared as a tool to construct free analogs of certain transport maps. Matrix-valued stochastic processes as solutions of a fSDE was the point of study in [16]. Evolution equations in non-commutative probability can be found in the dissertation [17]. For further works on fSDEs see [12], [13], [16], and [17].
Similar to the classical case, it is understandable that solutions of fSDEs may not be found explicitly. Hence, numerical methods have come into play to obtain approximation solutions to the underlying fSDE. At first, one can think of 1 as an equation in an abstract \(M_N(\mathbb{C})\)-valued algebra. But by the previously mentioned approximations of the distribution of elements \(X\in\mathcal{A}\) by matrix-valued elements \(X^{(N)}\), one can consider 1 in an abstract von Neumann algebra \(\mathcal{A}\). It is then natural to approximate the solutions of free stochastic differential equations by numerical methods on a matrix level. Looking at the classical counterpart, we have the Euler-Maruyama as well as the general Milstein scheme at hand ([18], [19]), where the latter acts as prototype for attempting first order methods. One of the major questions concerns the speed of convergence of the numerical iterations and here especially in the strong sense. In [20] the authors developed a free analog of the Euler-Maruyama scheme (fEMM) converging of order \(\gamma=1/2\) in the strong sense in \(L_1(\varphi)\) and order one in weak sense. These results were extended to defined stochastic theta methods for fSDEs and studying numerical stability ([21]).
In this paper, we first extend the results on the Euler-Mayurama scheme to general \(L_p(\varphi)\)-spaces (\(1\le p\le \infty\)). The reason is, that the free analog of the Burkholder-Gundy inequality is also valid in \(L_\infty(\varphi)\) ([3]). Second, we are seeking for Milstein type methods of order \(\gamma=1\) in the strong sense. Similar to the classical case the major ingredient in developing higher order methods is a stochastic Itô-Taylor expansion of the functions of drift and diffusion terms of the underlying fSDE. Here, the deep result of N. Azamov, A. Carey, P. Dodds and F. Sukochev ([22]) on multiple operator integrals (MOI) hits the scene. For details about MOIs, we also refer to the book [23]. The Itô-Taylor approximation and the representation of the derivatives in operator sense enables us to formulate a iterated stochastic Itô-Taylor expansion of the solution of the fSDE in order to identify those terms leading to methods of strong order \(\gamma=1\). We show that terms of higher order in the Itô-Taylor expansion are represented using multiple, iterated operator integrals. As in the commutative case, to gain speed of \(\gamma=1\) in the strong sense, such iterated integrals and the MOIs need to discretized. Nevertheless, in the non-commutative setting the situation is therefore more complex than for ordinary stochastic differential equations.
In the commutative case, when the SDE is equipped with a single Brownian motion \(B_t\in\mathbb{R}\), the double iterated integrals can be calculated directly as \(\int_0^t\int_0^sdB_u=\frac{1}{2}\left(B_t^2 -t\right)\) ([18]). For commutative SDEs with multiple Brownian motion, a direct calculation of the iterated integrals is possible for SDEs with special structure, so-called SDEs with commutative noise [18]. For other types of SDEs, the iterated integrals can only be approximated. We refer to [18], [24] and the references therein for an overview on approximation methods. In the non-commutative setting, in order to calculate iterated integrals directly, first, the Itô formula in product form ([3] has to be imposed and second, the MOIs need to be discretized. The Itô formula allows to resolve iterated free stochastic integrals into a product. We will see, that due to non-commutativity, the applicability of the product formula to the terms of the Itô-Taylor expansion is necessarily limited to the case of a single diffusion term, i.e. \(d=1\). But it turns out, that additionally, to be able to apply the free Itô formula, it is required to commute factors in the MOIs. The commuting operation must be handled by certain perturbation formulas (e.g. see [23]), but it comes with the costs of the appearance of triple operator integrals representing the second derivative of the diffusion term. These extra terms do not appear in the commutative setting.
Therefore a direct calculation of the iterated integrals is only possible in two cases. First, if the functions of the diffusion terms are affine functions of the unknown \(X_t\), i.e. the triple operator integrals vanish. Second, if one restricts to the case of \(\mathcal{A}=\mathcal{M}^N(\mathbb{R})_{sa}\), the triple operator integrals can be calculated directly with the help of the spectral distribution of \(X_t\) (resp. of the numerical approximation \(\overline{X}_t\).).
A numerical method for the general nonlinear case with multiple diffusion terms (\(d>1\)), requires the approximation of the free iterated integrals. We do this by resolving the integrals on a finer discretization of step size \(\delta t\) on each discretization interval \(\Delta t\) of the fSDE, where \(\delta t \leq \Delta t^2\). We develop the method (fSM) starting from the iterated Itô-Taylor expansion and discretize the terms necessary to achieve first order convergence. We will show, that the method can be simplified for \(d=1\). The approximation of the iterated integrals on the subintervals with length \(\delta t\) do not require the discretization of the previously mentioned triple operator integrals (which in fact, is unknown). Finally we give the proof of strong convergence order of \(\gamma=1\). The generality of the method comes with extra computational effort. All together, to carry over the Milstein scheme ([19], [18]) to the non-commutative case is only possible for \(d=1\) and affine diffusion coefficient. In this special case, it is possible to obtain a derivative free, first order method, where the iterated integrals are resolved exactly. For \(d=1\) and nonlinear diffusion, the derivatives of operator functions can be calculated from the spectral distribution of the solution in matrix level. This allows for a Milstein type method in the non-commutative setting, which does not appear in the classical commutative case. Finally we give numerical examples to reproduce the theoretical results numerically on matrix level both for \(d=1\) and \(d>1\).
The paper is organized as follows. 2 and 3 contains some preliminaries on free stochastic differential equations. 4.2 presents the technique on multiple operator integrals, which in fact serves for the Taylor-like expansion. We give some technical lemmas, which are required in the following. 5 is mainly a summary of [20] for the definition of strong convergence and the Euler-Maruyama method. It also states the first step of the Itô-Taylor expansion. 6 then shows the second iteration of the Itô-Taylor expansion and the necessary steps to develop a method of order \(\gamma=1\). In 19 we formulate a condition on the numerical method to be of order \(\gamma=1\). 20 gives a proof to the statement, that on finite time intervals and given any convergent numerical method in the strong sense (\(\gamma>0\)), the numerical is uniformly bounded on finite time intervals. 7 deals with the design of methods of order \(\gamma=1\). We first turns to those cases, for which the iterated integrals can be resolved directly. It is followed by the construction of methods for the general case \(d>1\). 8 is devoted to several examples to show the desired convergence rates numerically.

2 Preliminaries - Free Stochastic Calculus↩︎

Consider a classical probability space \((\Omega, \mathcal{F}, \mu)\) and random variables as measurable functions \(X:\Omega \rightarrow \mathbb{R}\). By taking an algebraic viewpoint, one can consider the algebra of random variables and their expectations as a fundamental concept. It allows to generalize the classical, commutative probability to more general settings, e.g. where the random variables are non-commutative. The space \(\mathcal{M}^N(\mathbb{C})=L_\infty\left(\Omega,\mu,\text{M}_N(\mathbb{C})\right)\) builds up a \(*\)-algebra with the unit matrix as identity and \(\varphi(M)=\frac{1}{N}\mathbb{E}(\text{tr}(M))\) as a trace mapping the identity matrix to \(1\). The classical notion of expectation is then replaced by the linear function \(\varphi\) on \(\mathcal{M}^N(\mathbb{C})\).
Free probability was created by C. Voiculescu in mid 1980’s by studying properties of von Neumann algebras. He introduced the notion of freeness, which extends the notion of independent random variables to non-commutative setting. He also discovered, that random matrices satisfy the freeness conditions asymptotically. By the help of non-commutative algebras it is possible to develop non-commutative probability theory. The limits \(N\rightarrow\infty\) can be handled properly in algebraic structures and lead to fruitful concepts. It turns out that non-commutative probability theory is realized by using operator algebras such as von Neumann algebras. We refer e.g. to [8], [25], [26] and references therein, for setting up non-commutative probability theory and relations to random matrices. To be complete, we give the following general definition (see e.g. [8]).

Since we consider von Neumann algebras with a unital, faithful and normal trace \(\varphi:\mathcal{A}\rightarrow\mathbb{C}\), we can introduce for \(1\leq p<\infty\) a norm on \(\mathcal{A}\) by \(\|X\|_p=\varphi(|X|^p)^{\frac{1}{p}}\). The Banach space completion of \(\mathcal{A}\) by \(\|\cdot\|_p\) is denoted by \(L_p(\varphi)\) (see e.g. [27]). Since the trace is finite we may consider \(\mathcal{A}\) as a subset of the predual \(L_1(\varphi)\) of the von Neumann algebra \(\mathcal{A}=L_\infty(\varphi)\). By \(\| \cdot \|\) we denote the usual operator norm in \(\mathcal{A}\). Let \(\mathcal{A}^{sa}=\{a\in\mathcal{A}, a^*=a\}\). For a non-commutative random variable \(X\in\mathcal{A}^{sa}\), there is a unique probability measure on \(\mathbb{R}\) with compact support having the same moments as \(X\) (e.g. [8], [26]). This probability measure is the distribution of the non-commutative random variable X.
The non-commutative analog of independence of classical random variables is the concept of free independence, or shortly freeness, of subalgebras of \(\mathcal{A}\) ([8]). Let \(\mathcal{A}_1,\dots \mathcal{A}_n\) be a family of \(n\in\mathbb{N}\) subalgebras of \(\mathcal{A}\). They are called freely independent (or simply free) in the sense of Voiculescu, if \(\varphi\left(X_1X_2\dots X_m\right)=0\) whenever the following conditions

1. \(X_j\in\mathcal{A}_{i(j)}\), where \(i(1)\neq i(2), i(2)\neq i(3), \dots , i(n-1)\neq i(n)\), \(j=1,\dots,m\)

2. \(\varphi(X_i)=0\) for all \(i=1,\dots,n\)

hold ([8]). If \(X\in\mathcal{A}\) is a self-adjoint element, then there is a unique spectral measure \(\mu\) on \(\mathbb{R}\) so that the moments of \(X\) are the same as the moments of the probability measure \(\mu\) defined by \(\varphi(X^k)=\int_\mathbb{R}x^k d\mu(x),\) see [8]. An important role plays the Cauchy transform \(G_X\) of \(\mu\) defined by \(G_X(z)=\int_{\mathbb{R}} \frac{d\mu(x)}{x-z},\) which is an analytic function defined on \(\mathbb{C}^+\) with values in \(\mathbb{C}^+\). The Cauchy transform \(G_X\) is the expectation of the resolvent of \(X\), i.e. \(G_X(z)=\varphi\left(\left(X-z\right)^{-1}\right)\).
The Cauchy transform carries all the properties of the spectral probability distribution of the self-adjoint operator \(X\). In [5] the authors used the Hilbert transform of the solution of the underlying fSDE to obtain detailed information about the distribution of the solution. In kargin? these results were extended to general fSDEs. They can be handled by a corresponding deterministic partial differential equations of the Cauchy transform \(G_X\). We will strongly depend on these results since it allows us to check the numerical results.

2.1 Free Brownian Motion↩︎

Motivated by the concept of classical Brownian motion the definition within non-commutative probability is as follows. Consider a von Neumann algebra \(\mathcal{A}\) with a faithful normal trace \(\varphi:\mathcal{A}\rightarrow \mathbb{C}\). A filtration \(\mathbb{F}=(\mathcal{A}_t)_{t \geq 0}\) is a family of von Neumann subalgebras \(\mathcal{A}_t\) of \(\mathcal{A}\) with \(\mathcal{A}_s \subset \mathcal{A}_t\) for \(s\leq t\). A family of elements \((X_t)_{t\geq 0}\subseteq \mathcal{A}\) is called adapted to the filtration \(\mathbb{F}\), if \(X_t\in\mathcal{A}_t\) for all \(t\geq 0\).

On finite time intervals \(J\subset \mathbb{R}\) the elements \(W_t\) are uniformly bounded, i.e. \(\sup_{t\in J}\|W_t\|<\infty\).

2.2 Stochastic Integration with Respect to Free Brownian Motion, Free Itô Formula↩︎

So called "free stochastic integrals", i.e. stochastic integrals with respect to free Brownian motion were introduced in [1], [3], [6]. Due to non-commutativity the stochastic integral is build up by integrands, where operators are multiplied both on the left and right of the integration variable. Such integrals are constructed by first introducing piecewise constant processes as integrands, so called simple biprocesses. By defining an appropriate norm, the vector space of simple biprocesses can be completed to the general space of biprocesses. We briefly summarize the construction of free integrals, for details we refer to [3].
Consider the opposite algebra \(\mathcal{A}^{op}\) to \(\mathcal{A}\) and a decomposition \(0=t_0<t_1<\dots t_m<\infty\). A simple adapted biprocess \(U_t\) is a piecewise constant map \[U_t= \begin{cases} A_k\otimes B_k,&t_k\leq t < t_{k+1}\\ 0,&t_n\leq t \end{cases},\] where \(A_k,B_k\in\mathcal{A}_{t_k}\) according to the filtration \(\mathbb{F}\). Following [3], the stochastic integral of \(U\) with respect to the Brownian motion \((W_t)_{t\geq 0}\) is the integral (\(\Delta W_k=W_{t_{k+1}}-W_{t_k}\)) \[\int_0^{\infty} U_s \sharp d W_s =\int_0^\infty A_sdW_sB_s := \sum_{k=0}^{m-1} U_{t_k} \sharp\left(\Delta W_k\right)=\sum_{k=0}^{m-1} A_{t_k}^j\left(\Delta W_k\right) B_{t_k}^j.\] By bilinearity, one can expand this definition to finite sums \(U_t=\sum_{i=1}^n A^i_t\otimes B^i_t\). An Itô-isometry can be obtained for all adapted simple biprocesses. E.g. for \(U=A\otimes B\) and \(V=C\otimes D\) the isometry reads \[\varphi\left[\int U_t \sharp d S_t \cdot\left(\int V_t \sharp d S_t\right)^*\right]=\int\left\langle U_t, V_t\right\rangle_{L_2(\varphi) \otimes L_2(\varphi)} dt = \int \varphi(A_tC_t^*)\varphi(B_tD_t^*)dt .\] The vector space of simple biprocesses can be equipped with the norms \[\|U\|_{\mathcal{B}_p}:=\left(\int\left\|U_t\right\|_{L^p\left(\varphi \otimes \varphi^{o p}\right)}^2 dt\right)^{1 / 2}.\] In the case of a selfadjoint biprocess \(U=A\otimes B\) the norm is \[\|A\otimes B\|_{\mathcal{B}_p}= \left(\int\varphi(A_t^2)\varphi(B_t^2)dt\right)^{1/2}.\] This vector space can be completed to the space \(\mathcal{B}_p\) (resp. for \(\mathcal{B}_p^a\) for the closed subspace of adapted processes). This allows to extend the mapping \(U\mapsto \int U_t\#dW_s\) isometrically to \(\mathcal{B}_2^a\rightarrow L_2( \varphi)\). Free calculus has the essential property, that free stochastic integrals are bounded operators even for \(p=\infty\). In the following, we proceed adapted to our needs. Let \(A_t, B_t\in\mathcal{A}\) and \((A_t)_{t\geq 0}, (B_t)_{t\geq 0}\) be continuous processes adapted to the filtration \(\mathbb{F}\). We use the simple notation of the free integrals as \(\int_{o}^{t} A_s dW_s B_s.\) As stated above, in non-commutative setting it is possible to set up a Burkholder-Gundy martingale inequality even for \(p=\infty\) ([3]): \[\label{ineq:BurkholdGundy} \left\| \int_0^t U_s\#dW_s\right\| \leq 2\sqrt{2}\left( \int_0^t \|A_s\|^2 \|B_s\|^2ds \right)^{\frac{1}{2}}.\tag{2}\] Hence, the inequality implies \(\left\| \int_t^{t+\Delta t} A_sdW_sB_s \right\|=O(\sqrt{\Delta t})\).

2.3 Free Itô-Process↩︎

Note, that the product \(b^i(s)dW_sc^i(s)\) is not necessarily self adjoint. We require \(X_t\in\mathcal{A}^{sa}\), so \(b^i, c^i\) cannot be chosen arbitrarily. It would be natural to define the diffusion terms symmetric, i.e. \(\sum_{i=1}^d b_t^idW_tb_t^i\). Considering the diffusion \(\int_0^tX_s dW_s + dW_sX_s\), it would be possible to write \(X_s dW_s + dW_sX_s = (X_s+\text{id})dW_s(X_s+\text{id}) + (X_s)(-dW_s)(X_s) - dW_s\). This requires mixed positive and negative signs before the free Brownian motion at different summands and would make the formalism more complicated. If one summand shows \(b^i\neq c^i\), there must be another summand \(c^idW_s b^i\) in order to meet the self adjoint requirement.
To keep it simple, we continue with the formalism ?? . For example, to express \(\int_0^tX_s dW_s + dW_sX_s\), we set \(d=2\) and \(b^1(X_t)=c^2(X_t)=X_t\) and \(b^2(X_t)=c^1(X_t)=\text{id}\).

3 Free Stochastic Differential Equations (fSDEs)↩︎

From above definition it follows immediately, that locally operator Lipschitz functions are continuous.

We take the following local existence and uniqueness results from kargin?.

As an initial example consider the free analog of the Ornstein-Uhlenbeck process (see kargin?) defined by the fSDE \[\label{ornsteinuhlenbeck} dX_t=\theta X_tdt + \sigma dW_t, \, t\geq 0, \,\,\theta, \sigma\in\mathbb{R}.\tag{3}\] Spectral information about the solution can be obtained by taking the Cauchy transform \(G\) of the self-adjoint element \(X_t\). \(G\) fulfills a deterministic partial differential equation (kargin?). Applying the Stieltjes inversion formula (see kargin?) to its solution, it is possible to recover the spectral distribution of \(X_t\). In the case \(\theta<0\) it turns out that the density of \(X_t\) is a semicircle distribution with radius \[R(t)=\sqrt{\frac{2\sigma^2}{|\theta|}(1-e^{-2|\theta| t})}.\] For \(t\rightarrow \infty\) the probability distribution function (PDF) converges to a semicircle with radius \(\sigma \sqrt{\frac{2}{|\theta|}}\). The case \(\theta\geq0\) is treated in the same way.

For more examples and discussion about the properties of the spectral distribution of the solution we refer to kargin?. Important fSDEs with a single diffusion term are

• Geometric Brownian Motion: \[dX_t=\theta X_tdt + \sqrt{X_t}dW_t\sqrt{X_t}\] By applying \(\varphi\) it follows \(\varphi(X_t)=e^{\theta X_t}\).

• As shown in kargin?, the solution of the following equation explodes in finite time \[dX_t=kX_tdW_tX_t.\] If \(X_0=I\), then the spectral distribution of \(X_t\) is defined for all \(t\leq 1/k^2\).

Equations with a single diffusion term

• Geometric Brownian Motion 2: \[dX_t=\theta X_tdt + X_tdW_t + dW_tX_t\] This is an example of as fSDE with two diffusion terms. The Brownian motion is the same in both summands.

• Free CIR Equation from [14], \[dX_t=(a-bX_t)dt + \sigma/2(\sqrt{X_t}dW_t + dW_t\sqrt{X_t}).\]

 
 

4 Operator Integrals, Taylor Series, Itô-Formula↩︎

In this section we introduce some well-known but necessary ingredients. We start with

4.1 Operator Integrals↩︎

An important ingredient in the following is the notion of a multiple operator integral, originating in the study of perturbations of operator functions. For functions with special properties, it is possible to express the difference \(f(X)-f(Y)\), where \(X,Y\in \mathcal{A}\), as a so-called double operator integral. This idea can be extended, such that under certain conditions, the Fréchet-derivatives of operator functions can be formulated as operator integrals. We refer to [23] and [22].
The definition of a double operator integrals \(T_{f^{[1]}}\), resp. a triple operator integral \(T_{f^{[2]}}\) is given in [22] and [22]. Recalling: For \(X\in\mathcal{A}^{sa},\,Y\in\mathcal{A}\), \[\label{def:optintgeneral} \begin{align} T_{f^{[1]}}^{X,X}(Y) &= \int_{\Pi^{(2)}}e^{i(s_0-s_1)X} Y e^{is_1X}d\nu_f^{(2)}(s_0,s_1), \\ T_{f^{[2]}}^{X,X,X}(Y,Y) &= \int_{\Pi^{(3)}} e^{i(s_0-s_1)X} Y e^{i(s_1-s_2)X} Y e^{is_2X} d\nu_f^{(3)}(s_0,s_1,s_2). \end{align}\tag{4}\] The definition of the set \(\Pi^{(n)}\) and the measure \(\nu_f^{(n)}\) can be found in [22].
In our case, we consider functions, taken from the so-called Wiener space \(W_n(\mathbb{R})=\{f\in C^n(\mathbb{R}):~f^{(k)}, \mathcal{F}f^{(k)}\in L_1(\mathbb{R}), \, k=0,\dots,n\}\), \(n\in\mathbb{N}\). If \(f\in W_1(\mathbb{R})\) the double operator integral is well defined (resp. \(f\in W_2(\mathbb{R})\) for the triple operator integrals.) As stated in [4], \(C^2\) functions are locally operator Lipschitz. Since \(W_n(\mathbb{R})\subset C^2(\mathbb{R})\) for \(n\geq 3\), functions \(f\in W_n(\mathbb{R}), n\geq 3\) are locally operator Lipschitz in operator norm. P. Biane and R. Speicher ([3]) developed free calculus with respect to free Brownian motion.
For functions \(f\in W_n(\mathbb{R})\), it is possible to give a Taylor approximation with appropriate remainder term [23], [22]. For example, let \(f \in W_3(\mathbb{R})\). Applying the Taylor series expansion [22], then the estimation of the remainder [23] yields \[\label{eq:uuu1} f(X_{i+1})-f(X_i) = T_{f^{[1]}}^{X_i,X_i}(\Delta X) + T_{f^{[2]}}^{X_i,X_i,X_i}(\Delta X,\Delta X) + O(\|\Delta X\|^3)\tag{5}\] where \(\Delta X=X_{i+1}-X_i\) and \(X_{i+1},X_i\in\mathcal{A}\).
The following lemmas are needed in the following. Although they are well known in the literature, we state these lemmas for better readability. [lem:referee] states that e.g. if a left factor \(X\in\mathcal{A}\) is pushed into the operator integral, then due to non-commutativity, a triple operator integral appears. 10 allows to replace a double operator integral by a difference of the underlying function, together with an estimation of the reminder. It allows to discretize the double operator integral to be implemented in a numerical method. 11 and 13 are needed e.g. in the proof of 20.

The following lemma is a simplified version of [23] resp. [22]. For our purposes, it is sufficient to consider functions \(f\in W_n(\mathbb{R})\). Since \(W_n(\mathbb{R})\subset B_{\infty 1}^n\) the assumptions of [23] are fulfilled (also for the case \(p=\infty\)). For the definition of the Besov space \(B_{\infty 1}^n\) we refer to [23].

Note that we rely mainly on the norm estimation of the remainder \(R_{2,f,A}(V)\) but not on the specific expression.

4.2 Free Itô Formula↩︎

An important ingredient in the development of numerical methods for fSDEs and their convergence properties is a free analog of the Itô formula. Influenced by the formalism from the commutative Wagner-Platen extension for solutions of SDEs ([18]), we adapt the notation of the Itô-formula introduced in [3].

For the proof we refer to [3]. In [9], the theorem was derived directly from a Taylor expansion of operator function \(f\).
In the context of stochastic differential equations, we will make use of 14 to work out an iterated Itô-Taylor expansion of the solution of the fSDE (see introduction of 5 and 7.2.1).
The Itô formula can also be expressed in product form, see [3], kargin?: \[\label{eq:ItoProductIntegral} \begin{align} \int_0^1a_tdW_tb_t\int_0^1c_tdW_tdt &=\int\left( \int_0^t a_sdW_sb_s\right)dW_tc_tdW_td_t\\ &+\int a_tdW_tb_T \left( \int_0^t a_sdW_sb_s\right)\\ &+\int_0^1 \varphi(b_tc_t)a_td_tdt. \end{align}\tag{7}\] This formula will play a major role in developing first order methods. The rule can also be written in differential form as (see kargin?) \[\label{freeItoFormulaAsInKargin} a_tdW_tb_t \cdot c_tdW_td_t= \varphi\left(b_tc_t\right)a_td_tdt.\tag{8}\] In the important case \(a_t=c_t=d_t=1\) this yields the formal rules \[dW_tb_tdW_t = \varphi\left(b_t\right)dt,\quad dW_tdW_t=dt.\] The formula also implies the formal rules \[dt^2=dW_tdt=0.\]

5 Numerical Approximation, Euler Maruyama↩︎

Let \(0<T<\infty\) and \(I=[0,T]\). We assume that the fSDE ?? has a unique solution \((X_t)_{t\in[0,T]}\) adapted to the filtration \(\mathbb{F}\). For \(L\in\mathbb{N}\), consider a discretization of \(I\) as \(t_0=0<t_1<\cdots<t_L=T\) and \(\Delta t=T/L\). We recurrently calculate a numerical approximations \((\overline{X}_k)_{k=0,\dots,L}\), \(\overline{X}_k\in\mathcal{A}_{t_k}\), by \[\overline{X}_{k+1}=\overline{X}_k + \Phi\left(t_k,\overline{X}_k, \Delta t, W_{t_{k+1}}-W_{t_k}\right),\] starting with \(\overline{X}_0=X_0\in\mathcal{A}^{sa}\). \(\Phi:\mathbb{R}\times \mathcal{A} \times \mathbb{R}\times \mathcal{A} \rightarrow \mathcal{A}\) is the increment function, which determines the numerical algorithm. In the following we simplify the notation, and write \(X_k\) as the solution of ?? evaluated at time point \(t_k\) i.e. \(X_k=X_{t_k}\). Analog for \(W_{t_k}=W_k\), \(\Delta W_k=W_{{k+1}}-W_{k}\).

In [20], the authors introduced the free analog of the Euler-Maruyama method, and proved strong convergence in \(L_2(\varphi)\) of order \(\gamma=1/2\) (together with weak convergence of order \(\gamma=1\)). In this paper, we extend this result to strong convergence in \(L_\infty(\varphi)\) with order \(\gamma=1/2\). This result will follow immediately from considerations towards methods of order \(\gamma=1\). Nevertheless, we repeat the definition of the free Euler-Maruyama (fEMM) in the following for the sake of completeness. In 5.1, we perform the first step of the stochastic Itô-Taylor expansion, already carried out in [20], but adding 16. In 5.2 the method fEMM will be defined together with the statement of the convergence result in \(L_\infty(\varphi)\).

5.1 Stochastic Itô-Taylor Expansion - First Step↩︎

Assuming \(a,b^i,c^i\in W_3(\mathbb{R})\) in ?? , we can apply the free Itô formula ?? for \(f=a,b^i,c^i\). This yields an iterated free Itô formula which allows to motivate and define a free analog of the Euler-Maruyama method (see also [20]). Using the abbreviations \(a(X_t)=a_t\) (similar notation for \(b^i,c^i\)), we obtain \[\begin{gather} \label{eq:h1231} X_{t+\Delta t}-X_t=\int_t^{t+\Delta t} a_t ds+ \int_t^{t+\Delta t}\int_t^s L^0[a_u]du\,ds+\int_t^{t+\Delta t}L^1[a_u]_t^{s}ds + \\ \sum\limits_{i=1}^d\int\limits_t^{t+\Delta t}\left( b_t^i+\int_t^sL^0[b_u^i]du+L^1[b_u^i]_t^{s}) \right)dW_s\left( c_t^i+\int_t^sL^0[c_u^i]du+L^1[c_u^i]_t^{s}) \right) \end{gather}\tag{9}\] Since \(a_t,b_t^i,c_t^i\) do not depend on the integration variable \(s\), we rewrite 9 as \[\label{eq:h4} X_{t+\Delta t}-X_t=a_t\Delta t+\sum_{i=1}^db_t^i(W_{t+\Delta t}-W_t)c_t^i + \sum_{i=1}^d M^i_t(\Delta t) + \sum_{i=1}^dR_t^i(\Delta t),\tag{10}\] where \[\label{eq:h5-1} M_t^i(\Delta t) = \int_t^{t+\Delta t} b_t^i dW_s\left(L^1[c_u^i]_t^{s} \right)+ \int_t^{t+\Delta t}\left(L^1[b_u^i]_t^{s}\right)dW_sc_t^i\tag{11}\] and \[\begin{gather} \label{eq:h5} R_t^i(\Delta t) = \int_t^{t+\Delta t}\int_t^s L^0[a_u^i]du\,ds+ \int_t^{t+\Delta t}L^1[a_u^i]_t^sds + \\ +\int_t^{t+\Delta t} b_t^i dW_s\left( \int_t^sL^0[c_u^i]du\right) + \int_t^{t+\Delta t}\left(\int_t^sL^0[b_u^i]du\right) dW_s c_t^i + \\ +\int_t^{t+\Delta t}\left(\int_t^sL^0[b_u^i]du\right) dW_s\left( \int_t^sL^0[c_u^i]du\right)+\\ +\int\limits_t^{t+\Delta t}\left(\int_t^sL^0[b_u^i]du\right) dW_s\left( L^1[c_u^i]_t^{s} \right) +\int\limits_t^{t+\Delta t}\left(L^1[b_u^i]_t^{s}\right)dW_s\left( \int_t^sL^0[c_u^i]du\right) + \\ +\int_t^{t+\Delta t}\left(L^1[b_u^i]_t^{s}\right)dW_s\left( L^1[c_u^i]_t^{s} \right). \end{gather}\tag{12}\]

Due to the boundedness and continuity of the involved functions \(a,b^i,c^i\) the above integrals are well defined.

5.2 Strong Convergence of the free Euler-Maruyama Method↩︎

The free Euler-Maruyama method can now be motivated from 10 simply by skipping the terms \(M_t^i\) and \(R_t^i\). The following is taken from [20].

In [20] the authors proved strong convergence order of \(\gamma=\frac{1}{2}\) of fEMM ?? in \(L_2\)-norm under certain assumptions. In [21] the authors developed implicit methods supplemented by stability analysis.
We are now ready to formulate the following theorem, which extends the \(L_2(\varphi)\) convergence of fEMM to \(L_\infty(\varphi)\) in the strong sense. For better readability, we formulate the theorem already at this point, although it is a direct consequence of 19. We refer the reader to 6 for details.

We refer to [cor:thm1-FEMM-order-1] which states strong convergence of fEMM of order \(\gamma=1\) in case the diffusion is constant, i.e. for equations of the form \(dX_t=a(X_t)dt+ \lambda dW_t, \lambda \in \mathbb{R}\).

6 Approximations of Strong Order \(\gamma=1\) in \(L_p(\varphi), 1\leq p \leq \infty\)↩︎

Starting point is the expansion 10 , which includes 12 and 11 . By iterating the expansion (in 6.1), we will obtain those terms, which must be discretized in order to obtain a method with \(\gamma=1\). We will see, that these terms contain iterated free stochastic integrals, just as in the commutative case. A proper discretization of the iterated integrals is required, which will be content of 7. Unfortunately, it turns out, that only in very special cases, the iterated integrals can be calculated directly. Only the cases

• single affine diffusion term (\(d=1\)) in general \(\mathcal{A}\),

• single nonlinear diffusion term in the special von Neumann algebra \(\mathcal{A}=\mathcal{M}^N(\mathbb{R})_{sa}\),

end up in a derivative free Milstein-type method for fSDEs. All the other cases do rely on a approximation of the iterated integrals. In the case of multiple nonlinear diffusion terms, the iterated integrals in \(m_t(\Delta t)\) cannot be resolved by Itô (see 7.1.3). Up to now, the general case \(d>1\) can only be handled by approximating the iterated integrals via the method of subdividing the discretization interval.

6.1 Stochastic Itô-Taylor Expansion - Second Iteration↩︎

As in the commutative case, the key in the development of higher order methods is the iterated stochastic Itô-Taylor expansion of the solution of the underlying differential equation and its discretization. To obtain such an expansion, one needs to apply the free Itô formula ?? repeatedly to the drift and diffusion terms \(a, b^i, c^i\) of the underlying fSDE. In the last chapter we defined fEMM out of the expansion 10 . To be able to define a method of strong order \(\gamma=1\), we apply ?? to \(a,b^i,c^i\) in 11 once more, but we only take the terms \(b_t^i\) and \(c_t^i\), which are the functions \(b^i,c^i\) evaluated at time point \(t\), for example \(b^i(X_t)=b_t^i\). This yields \[\label{eq:proof-strong-conv-start-def} X_{t+\Delta t}=X_t + a_t\Delta t + \sum\limits_{i=1}^db_t^i\Delta W_t c_t^i + \sum\limits_{i=1}^dm_t^i(\Delta t)+ \sum\limits_{i=1}^d\rho_t^i(\Delta t),\tag{13}\] where the term \(m_t^i(\Delta t)\) contains the iterated free stochastic integrals, i.e. \[\begin{gather} \label{E:m1rep1-1-def} m_t^i(\Delta t) = \int_{t}^{t+\Delta t} b_t^i dW_s\left( T_{c^{i,[1]}}^{X_t,X_t}\left(\sum\limits_{j=1}^{d}b_t^j \int_{t}^s dW_uc_t^j \right)\right) + \\ + \int_{t}^{t+\Delta t}\left(T_{b^{i,[1]}}^{X_t,X_t}\left(\sum\limits_{j=1}^{d}b_t^j\int_{t}^sdW_uc_t^j\right)\right)dW_s c_t^i. \end{gather}\tag{14}\] By the same arguments as in 16, we have \(\|m_t^i(\Delta t)\|^2=\mathcal{O}(\Delta t^2)\). Note that \(m_t^i(\Delta t)\) are the only terms of \(\mathcal{O}(\Delta t^2)\) in the expansion 13 . All higher order terms (including \(R^i\) in 10 ) are collected in \(\rho_t^i(\Delta t)\). It follows \(\|\rho_t^i(\Delta t)\|^2 = \mathcal{O}(\Delta t^3)\) by the same arguments as for \(m_t^i\) above.

6.2 General Theorem of Strong Convergence with Order \(\gamma=1\) in \(L_p, \,1\leq p\leq \infty\)↩︎

Let \(t_0=0<t_1<\dots<t_L=T\) a discretization of the interval \([0,T]\). A numerical method based on the Itô-Taylor expansion is build up by finding a proper discretization of the terms in 14 , denoted by \(\overline{m}_k^i(\Delta t)\). The construction of such approximations is subject of 7. A numerical approximation \(\overline{X}_k\) to \(X_{t_k}\) is then calculated by \[\label{E:Milstein-Method-General} \overline{X}_{k+1}=\overline{X}_k+a(\overline{X}_k)\Delta t + \sum_{i=1}^d b^i(\overline{X}_k)\Delta W_kc^i(\overline{X}_k) + \sum_{i=1}^{d}\overline{m}_k^i(\Delta t).\tag{15}\] We will now formulate a general theorem, which imposes conditions on \(\overline{m}_k^i(\Delta t)\) in order to obtain a convergent method in the strong sense with order \(\gamma=1\).

In the special case of constant diffusion the term 14 vanishes and it follows

6.2.1 Boundedness of the Numerical Solution↩︎

In 19 we assumed a uniform bound on the numerical solution. We now show, that if a numerical method performs with strong convergence of order \(\gamma>0\), this assumption is true.

7 Construction of Methods↩︎

The development of a Milstein type method of order \(\gamma=1\) requires the discretization of the iterated free stochastic integrals 14 . It is natural to ask, under which circumstances these integrals can be calculated directly. As it turns out, there’s a major difference between \(d=1\) and \(d>1\). In the first case, Itô can be applied, in the second case not. Therefore this section starts with a discussion of the case \(d=1\). Due to the imposed non-commutativity the treatment of 14 differs from the commutative case even for a single diffusion \(d=1\). We have to deal with additional triple operator integrals, which stand in place for the second derivative of the diffusion operator functions. Unfortunately theses terms are only of order \(\Delta t^2\) (in square of \(L_\infty(\varphi))\)-norm, which does not allow to skip them in the development of first order methods. These terms do not appear in the commutative case, where the classical Milstein method only contains first order derivatives.
In 7.1 we discuss, under which conditions these operator integrals can be calculated directly in case of a single diffusion term. Unfortunately, we will see in section 7.2 that all these methods fail for the case \(d>1\). Even more, in 7.2 it is stated, that the Itô formula cannot be applied to resolve the iterated integrals as in the case \(d=1\). We therefore present a method to discretize the iterated integrals to achieve order \(\gamma=1\). In fact, to the best of our knowledge, there’s no other discretization of iterated integrals in the given context known in the literature.

7.1 Methods for Single Diffusion, \(d=1\)↩︎

Now we turn to the construction of methods 15 in the simplest case of a single diffusion term \(d=1\). We will see, that due to the non-commutativity, only in special cases a simple discretization can be found, as already mentioned.
In the commutative case, the key to the construction of higher order numerical approximation relies on the discretization of the iterated stochastic integrals in the stochastic Taylor expansion of the solution of the underlying SDE (e.g. see [18], [19]). Consider the simplest case of a scalar SDE with a single Brownian motion and differentiable diffusion term \(b\), i.e. \[dx_t=a(t,x_t)dt + b(t,x_t)dW_t.\] The stochastic Taylor expansion for a method of strong order \(1\) includes the following iterated integral, which can be integrated exactly, i.e. \[\int_{t_k}^{t_{k+1}}\int_{t_k}^s dW_u dW_s =\frac{1}{2}\left((\Delta W_{k})^2-\Delta t\right).\] The resulting method is the well-known Milstein method ([19], [18]). For SDE with more complicated structure, e.g. multiple Brownian motions, an efficient evaluation of the iterated integrals is only possible, if the structure of the SDE has special properties, i.e. commutative noise.
In the non-commutative case, the structure of the terms in the stochastic Taylor expansion 14 consists of a sum of iterated integrals. The goal is, to combine these summands in order to convert them into a product by applying the Itô formula in product form, [3].
It can easily be seen in 14 that the application of the product form of the Itô formula requires \(c^i=b^j, i,j=1,\dots d\). This means \(d=1\), which implies \(b=c\). In this case 14 simplifies to \[\begin{gather} \label{E:fMMvariante1-def-1} m_t(\Delta t)= b_t\left(\int_{t}^{t+\Delta t} dW_sT_{b^{[1]}}^{X_t,X_t}( \int_{t}^s b_tdW_u) \right. +\\+ \left. \int_{t}^{t+\Delta t} T_{b^{[1]}}^{X_t,X_t}(\int_{t}^s dW_u b_t )dW_s \right)b_t. \end{gather}\tag{24}\] Note that we applied the fact, that two elements \(U,V\in\mathcal{A}\) commute, if both belong to the same subalgebra, e.g. \(U,V\in \langle X_t\rangle \subset \mathcal{A}\).
In order to utilize the Itô formula we need to commute \(dW_s\) into the operator integrals, such that the outer stochastic integrals are pushed into the argument of the double operator integral \(T_{b^{[1]}}^{X_t,X_t}(\cdot)\). In the non-commutative case, this can be achieved by [lem:referee], but with extra costs. Applying [lem:referee] to 24 yields \[\begin{gather} \label{E:fMMvariante1-referee} m_t(\Delta t)= b_t\left(T_{b^{[1]}}^{X_t,X_t}( \int_{t}^{t+\Delta t} dW_s\int_{t}^s b_tdW_u) + \right. \\+ \left. T_{b^{[1]}}^{X_t,X_t}(\int_{t}^{t+\Delta t}\int_{t}^s dW_u b_tdW_s ) \right)b_t +\\+ b_t\left(T_{b^{[2]}}^{X_t,X_t,X_t}(V_t,U_t^l) - T_{b^{[2]}}^{X_t,X_t,X_t}(U_t^r,V_t)\right)b_t, \end{gather}\tag{25}\] where we shortened the notation by \(V_t=\int_t^{t+\Delta t}dW_sX_t-X_t\int_t^{t+\Delta t}dW_s\), and \(U_t^l=\int_t^{t+\Delta t} b_t dW_u\), \(U_t^r=\int_t^{t+\Delta t}dW_ub_t\). The first two summands in 25 are ready to be handled by Itô, i.e. \[\begin{gather} \label{E:fMMvariante1-def-2x} m_t(\Delta t)= b_t\left(T_{b^{[1]}}^{X_t,X_t}\left( \int_{t}^{t+\Delta t} dW_s b_t \int_t^{t+\Delta t}dW_s-\int_t^{t+\Delta t}\varphi(b_t)\boldsymbol{1}dt\right)\right)b_t +\\+ b_t\left(T_{b^{[2]}}^{X_t,X_t,X_t}(V_t,U_t^l) - T_{b^{[2]}}^{X_t,X_t,X_t}(U_t^r,V_t)\right)b_t \end{gather}\tag{26}\] Due to the additional triple operator integrals \(T_{b^{[2]}}^{X_t,X_t,X_t}\) we face new difficulties. Unfortunately due to 16, the two extra terms obey the property \[\|T_{b^{[2]}}^{X_t,X_t,X_t}(V_t,U_t^l)\|^2=\mathcal{O}(\Delta t^2),\,\|T_{b^{[2]}}^{X_t,X_t,X_t}(U_t^r,V_t)\|^2=\mathcal{O}(\Delta t^2).\] This shows, that these terms cannot be dropped. They request for a discretization in order to achieve a method of strong order \(\gamma=1\).

7.1.1 Linear Diffusion \(b\), \(d=1\)↩︎

Assume that the diffusion in the fSDE is linear, i.e. \(b(X_t)= X_t\) (a factor \(\mu\neq 0\) will be imposed in the fSDE directly, see 22). This allows a direct calculation of the operator integrals in 26 . If \(b\) is affine, then \(T_{b^{[1]}}^{X,X}(Y)=Y , \, X,Y\in\mathcal{A}\). Since the divided difference \(b^{[2]}=0\), the terms including \(T_{b^{[2]}}^{X_t,X_t,X_t}\) in 26 vanish. We then obtain \[m_t(\Delta t)=X_t\left(\Delta W_t X_t \Delta W_t - \varphi(X_t)\Delta t\right) X_t.\] The discretization is the simply constructed via \[\label{eq:deq1lindiffm} \overline{m}_k(\Delta t)= \overline{X}_k\left(\Delta W_k \overline{X}_k \Delta W_k - \varphi(\overline{X}_k)\Delta t\right) \overline{X}_t.\tag{27}\] In detail:

7.1.2 Nonlinear Diffusion - \(\mathcal{A}=\mathcal{M}^N(\mathbb{R})\), \(d=1\)↩︎

For a general von Neumann algebra \(\mathcal{A}\) and nonlinear \(b\) the discretization of \(m_t(\Delta t)\) in 26 is unknown, to the best of our knowledge. But if one restricts the considerations to the algebra \(\mathcal{A}=\mathcal{M}^N(\mathbb{R})_{sa}\), then it is possible to calculate the operator integrals by spectral analysis of \(X_t\) resp. \(\overline{X}_t\) directly. Therefore we can build up a Milstein Method of order \(\gamma=1\) based on the stochastic Taylor expansion.
In the non commutative probability space \(\mathcal{A}=\mathcal{M}^N(\mathbb{R})_{sa}\), the operator integrals can be expressed as ([9], [22], [23]) \[\begin{align} \label{E:Discret-OpInt-Matrix} T_{b^{[1]}}^{X,X}(Y)&=\sum_{\boldsymbol{\lambda} \in \sigma(X)^2} b^{[1]}(\boldsymbol{\lambda})P_{\lambda_1}YP_{\lambda_2},\\ T_{b^{[2]}}^{X,X,X}(Y_1,Y_2) &= \sum_{\boldsymbol{\lambda} \in \sigma(X)^3} b^{[2]}(\boldsymbol{\lambda})P_{\lambda_1}Y_1P_{\lambda_2}Y_2P_{\lambda_3}, \end{align}\tag{28}\] where \(\sigma(X)\) denotes the spectrum of \(x\in\mathcal{A}\) and \(P_{\lambda_i}\) are the projection matrices onto the eigenspace of the eigenvalue \(\lambda_i\in\mathbb{R}\). For readability we write 26 slightly shorter. Evaluated at \(t=t_k\) expression 26 reads \[\label{E:fMMvariante1-def-2x-1} m_k(\Delta t)= b_k\left(T_{b^{[1]}}^{X_k,X_k}\left( Y_k\right) + T_{b^{[2]}}^{X_k,X_k,X_k}(V_k,U_t^l) - T_{b^{[2]}}^{X_k,X_k,X_k}(U_k^r,V_k)\right)b_k,\tag{29}\] where \(Y_k=\int_{t_k}^{t_k+\Delta t} dW_s b_k \int_{t_k}^{t_k+\Delta t}dW_s-\int_{t_k}^{t_k+\Delta t}\varphi(b_{k})ds\). The numerical method can be constructed from 29 via \[\label{E:m-matrix-method-deq1} \overline{m}_k(\Delta t)= \overline{b}_k \left( T_{b^{[1]}}^{\overline{X}_k,\overline{X}_k}(\overline{Y}_k) +T_{b^{[2]}}^{\overline{X}_k,\overline{X}_k,\overline{X}_k}(\overline{V}_k,\overline{U}_k^l) - T_{b^{[2]}}^{\overline{X}_k,\overline{X}_k,\overline{X}_k}(\overline{U}_k^r, \overline{V}_k) \right)\overline{b}_k,\tag{30}\] where \[\begin{align} \overline{b}_k&=b(\overline{X}_{t_k})\\ \overline{Y}_k&=\Delta W_k \overline{b}_k \Delta W_k - \varphi(\overline{b}_k)\Delta t \\ \overline{V}_k&=\Delta W_k \overline{X}_k - \overline{X}_k \Delta W_k \\ \overline{U}_k^l &= \overline{b}_k\Delta W_k \\ \overline{U}_k^r &= \Delta W_k \overline{b}_k \\ \Delta W_k &= W_{k+1}-W_k. \end{align}\]

7.1.3 Nonlinear Diffusion - Subdivision Method, \(d=1\)↩︎

Now we drop the assumption \(\mathcal{A}=\mathcal{M}^N(\mathbb{R})_{sa}\) of 7.1.2 and consider the fSDE in a general von Neumann algebra \(\mathcal{A}\). Here, the problem is that for a method of strong order \(\gamma=1\), each of the terms including \(T_{b^{[2]}}^{X_t,X_t,X_t}\) in 26 have to be discretized. So far, there is no proper discretization of triple operator integrals \(T_{b^{[2]}}^{X_t,X_t,X_t}\) in a general von Neumann Algebra \(\mathcal{A}\) known (comparable to 10).
At this point it turns out, that in the non-commutative case, a literal free analog of the classical Milstein Method for commutative SDE’s with exactly evaluated iterated integrals seems to be impossible. Therefore in case of nonlinear diffusion \(b\) and a general von Neumann algebra \(\mathcal{A}\) it remains to apply the so called subdivision method, which we develop in the following. In 7.2 this method will be applied to the general case \(d>1\).

In the following we use the abbreviation \(\Delta W_{t,\tau}=W_\tau - W_t, \, 0\leq t\leq \tau\). If \(t\) and \(\tau\) are discretization points, then we only write the index of the discretization points as indices of \(\Delta W\) (same for \(X_t\)). Consider a subdivision of \([t,t+\Delta t]\) as described in 26. We use \(\Delta W_{l-1,l}=W_{\tau_{l}}-W_{\tau_{l-1}}\) frequently. If one of the discretization points coincide with the end points of the interval (\(\tau_0=t\)), we write the point in the index, e.g. \(\Delta W_{t,l-1}=W_{\tau_{l-1}}-W_t\).

7.2 Nonlinear Diffusion, \(d>1\)↩︎

The key in the development of the methods in 22 and 23 was the application of the Itô formula to resolve the iterated integrals to a product to obtain 26 . In the general case, we have to deal with \[\begin{gather} m_t^i(\Delta t) = \int_{t}^{t+\Delta t} b_t^i dW_s\left( T_{c^{i,[1]}}^{X_t,X_t}\left(\sum\limits_{j=1}^{d}b_t^j \int_{t}^s dW_uc_t^j \right)\right) + \\ + \int_{t}^{t+\Delta t}\left(T_{b^{i,[1]}}^{X_t,X_t}\left(\sum\limits_{j=1}^{d}b_t^j\int_{t}^sdW_uc_t^j\right)\right)dW_s c_t^i. \end{gather}\] To simplify the iterated integrals by the Itô formula in product form, it is algebraically required that \(c^i=b^j, i,j=1,\dots d\), but this is equivalent to \(d=1\). Unfortunately, the application of Itô is not possible for \(d>1\).
The general case \(d>1\) requires therefore the treatment of the iterated integrals. So far, to the best knowledge of the authors, there is no approximation method of non-commutative iterated free stochastic integrals known. Therefore we treat the integrals by the subdivision method developed in the previous chapter.

7.2.1 Subdivision Method, \(d>1\)↩︎

In this section we define a method of strong order \(\gamma=1\) for general \(\mathcal{A}\) and \(d>1\). It is based on approximating the iterated integrals in 14 by a proper subdivision of \(\Delta t\), such that, loosely speaking, the iterated integrals are approximated good enough to obtain a higher convergence rate of the numerical method.
As a first step, we simplify 14 . Since elements \(b^i(X_t), c^i(X_t)\) and \(e^{is X_t}\) belong to the same subalgebra generated by \(X_t\in\mathcal{A}^{sa}\), they do commute and due to the linearity of the operator integrals. We can rearrange 14 to \[\begin{gather} \label{E:m1rep1-1-def-1} \sum\limits_{i=1}^d m_t^i(\Delta t) = \sum\limits_{i=1}^d \sum\limits_{j=1}^d b_t^i\int_t^{t+\Delta t}dW_s T_{c^{i,[1]}}^{X_t,X_t}(b_t^j\int_t^sdW_u)c_t^j +\\+ \sum\limits_{i=1}^d\sum\limits_{j=1}^d b_t^j \int_t^{t+\Delta t}T_{b^{i,[1]}}^{X_t,X_t}(\int_t^sdW_u c_t^j)dW_sc_t^i =\\= \sum\limits_{i=1}^d \sum\limits_{j=1}^d b_t^i\left( \int_t^{t+\Delta t}dW_sT_{c^{i,[1]}}^{X_t,X_t}(b_t^j\int_t^s dW_u) \right. +\\+ \left. \int_t^{t+\Delta t}T_{b^{j,[1]}}^{X_t,X_t}(\int_t^sdW_uc_t^i) dW_s\right)c_t^j =\\= \sum\limits_{i=1}^d \sum\limits_{j=1}^d b_t^i\left(I_1^{i,j}(\Delta t) + I_2^{i,j}(\Delta t)\right)c_t^j \end{gather}\tag{33}\] Now consider a subdivision of \([t,t+\Delta t]\) as described in 26. We rewrite \(I_1^{i,j}\) and \(I_2^{i,j}\) using the step-width \(\delta t \leq \Delta t^2\), \[\begin{align} I_1^{i,j}(\Delta t) &= \sum_{l=1}^n\int_{\tau_{l-1}}^{\tau_l}dW_sT_{c^{i,[1]}}^{X_t,X_t}(b_t^j\int_t^s dW_u)\\ &= \sum_{l=1}^n\int_{\tau_{l-1}}^{\tau_l}dW_sT_{c^{i,[1]}}^{X_t,X_t}(b_t^j\int_{\tau_{l-1}}^s dW_u + b_t^j\int_t^{\tau_{l-1}}dW_u)\\ &=\sum_{l=1}^n\left[\int_{\tau_{l-1}}^{\tau_l}dW_sT_{c^{i,[1]}}^{X_t,X_t}(b_t^j\int_{\tau_{l-1}}^s dW_u)+\Delta W_{l-1,l}T_{c^{i,[1]}}^{X_t,X_t}(b_t^j\Delta W_{t,l-1})\right] \end{align}\]

Due to freeness we have the following estimation by applying [23] and considering that \(\left\|\Delta W_{t,t+\delta t} \right\|^2=\mathcal{O}(\delta t^2)\), \[\begin{gather} \label{eq:L2estimreferee} \|\int_{t}^{t+\delta t}T^{X_t,X_t,X_t}_{f^{[2]}}(dW_sX_t-X_tdW_s,b_t\int_t^sdW_u)\|^2 =\\=\|\int_{t}^{t+\delta t}T^{X_t,X_t,X_t}_{f^{[2]}}(\int_t^sdW_ub_t,X_tdW_s-dW_sX_t)\|^2 =O(\delta t^2). \end{gather}\tag{34}\] Now we are ready to apply [lem:referee] to the first summand of \(I_1^{i,j}(\Delta t)\). By [lem:referee] we push the integration variable \(dW_s\) into the operator integral \(T_{c^{i,[1]}}^{X_t,X_t}(\cdot)\) in the first summand. Together with 34 we obtain \[\begin{gather} I_1^{i,j}(\Delta t) =\sum_{l=1}^n\left[T_{c^{i,[1]}}^{X_t,X_t}(\int_{\tau_{l-1}}^{\tau_l}dW_sb_t^j\int_{\tau_{l-1}}^s dW_u) \right.+\\+\left.\Delta W_{l-1,l}T_{c^{i,[1]}}^{X_t,X_t}(b_t^j\Delta W_{t,l-1})\right] +A_1^{i,j}(\delta t). \end{gather}\] The term \(A_1^{i,j}(\delta t)\) is the sum over \(l=1,\dots,n\) of the triple operator integral produced through the application of [lem:referee]. Due to 34 it follows \(\|A_1^{i,j}(\delta t)\|^2 = \mathcal{O}(n\delta t^2)=\mathcal{O}(\Delta t^3)\), where the constant is independent of the index \(l\). Similarly we obtain \[\begin{gather} I_2^{i,j}(\Delta t) = \sum_{l=1}^n\left[T_{b^{j,[1]}}^{X_t,X_t}(\int_{\tau_{l-1}}^{\tau_l}\int_{\tau_{l-1}}^s dW_uc_t^idW_s) \right. +\\+\left. T_{b^{j,[1]}}^{X_t,X_t}(\Delta W_{t,l-1}c_t^i)\Delta W_{l-1,l}\right]+A_2^{i,j}(\delta t). \end{gather}\] As for \(A_1^{i,j}(\delta t)\) we conclude \(\|A_2^{i,j}(\delta t)\|^2=\mathcal{O}(\Delta t^3)\). Finally 33 is rewritten as \[\label{eq:mTermStochTaylor} \begin{align} m_t^i(\Delta t) &= \sum\limits_{j=1}^d b_t^i\left(I_1^{i,j}(\Delta t) + I_2^{i,j}(\Delta t)\right)c_t^j =\\ &= \sum\limits_{j=1}^d \sum\limits_{l=1}^n b_t^i\left[ T_{c^{i,[1]}}^{X_t,X_t}(\int_{\tau_{l-1}}^{\tau_{l}}dW_sb_t^j\int_{\tau_{l-1}}^s dW_u)\right]c_t^j + \\ &+ \sum\limits_{j=1}^d \sum\limits_{l=1}^n b_t^i\left[T_{b^{j,[1]}}^{X_t,X_t}(\int_{\tau_{l-1}}^{\tau_l}\int_{\tau_{l-1}}^s dW_uc_t^idW_s) \right]c_t^j + \\ &+\sum\limits_{j=1}^d \sum\limits_{l=1}^n b_t^i\left[ \Delta W_{l-1,l}T_{c^{i,[1]}}^{X_t,X_t}(b_t^j\Delta W_{t,l-1}) \right]c_t^j +\\ &+ \sum\limits_{j=1}^d \sum\limits_{l=1}^n b_t^i\left[T_{b^{j,[1]}}^{X_t,X_t}(\Delta W_{t,l-1}c_t^i)\Delta W_{l-1,l} \right]c_t^j\\ &+ \sum\limits_{j=1}^d (A_1^{i,j}(\delta t)+A_2^{i,j}(\delta t)). \end{align}\tag{35}\] In the general case \(d>1\) we have to deal with the condition \(b_t^j\neq c_t^i\). The consequence is, that the sum of the operator integrals \(T_{c^{i,[1]}}^{X_t,X_t}(\cdot)\) and \(T_{b^{j,[1]}}^{X_t,X_t}(\cdot)\) in 35 cannot be combined. Only the case \(d=1\) allows this sum to be simplified. Then, a simplification of iterated stochastic integrals by the Itô formula is possible. For details, we refer to 7.1 but for now, our aim is to define a numerical method suitable for the case \(d>1\). To do so, we now discretize 35 by 10 and skip the terms \(A_1^{i,j}\) and \(A_2^{i,j}\) in 35 .

8 Numerical Examples↩︎

In the following examples we present different fSDEs to confirm the insights from the previous chapters numerically. We select three fSDEs with \(d=1\), \(d=2\) and different functions \(b,c\). Additionally we show that an a posteriori estimation of the convergence order gives the expected results. We mainly follow [20] and numerically seek an approximation of the solution of the fSDE at a final time \(T>0\).

8.1 A Simple Example↩︎

We consider the fSDE \[\label{ex:explosive-sde} dX_t= X_tdW_tX_t\tag{42}\] with start value \(X_0=I\). In kargin? it is shown, that the spectral distribution of the solution \(X_t\in\mathcal{A}^{sa}\) exists of all \(t\leq1\) and is supported on the interval \[\left[\frac{(1-\sqrt{t})^2}{(1-t)^2},\frac{(1+\sqrt{t})^2}{(1-t)^2}\right].\] For \(t\in ]0,1]\), the density of the spectral distribution is given by \[\label{kap8-example-PDF} f(x) = \frac{\sqrt{-(1-T)^2 x^2+2(1+T)x-1}}{2\pi T x^3}.\tag{43}\] For \(t=1\) the density is supported on \([1/4,\infty[\).
We implement the method ?? with \(\lambda=1, \mu=0\) defined in 22. The probability density function (PDF) of the eigenvalues of the numerical solution are then an approximation of 43 .

The numerical algorithms fEMM and ?? are calculated with four different time steps \(\Delta t_R=\frac{T}{2^{13}}2^{R}\) with \(R=0,8,9,10\) and final time \(T=0.1\). The path of the free Brownian motion in the case \(R=0\) will be used to calculate all necessary moments for \(R=8,9,10\). The exact solution, which is unkown, will be simulated by \(R=0\). An a posteriori error estimation is shown in the succeeding examples. Let’s consider \(U\in\mathbb{N}\) different paths \(W_t(\omega_u), \, u=1,\dots,U\) of the free Brownian motion. In the following \(\overline{X}_T(u,\Delta t_R)\) denotes the numerical solution calculated at time point \(T\), by time step \(\Delta t_R\) and path \(W_t(\omega_u)\). We simply call \(W_t(\omega_u)\) as path number \(u\). For each time step \(\Delta t_R\) for \(R=8,9,10\) and each path number \(u\) the strong path-wise error in \(L_\infty(\varphi)\)-norm is calculated as \(e^{\text{Sim}}(u,\Delta t_R)=\|\overline{X}_T(u,\Delta t_R) - \overline{X}_T(u,\Delta t_0)\|\), where the numerical solution \(\overline{X}_T(u,\Delta t_R)\) is either calculated by fEMM or the method ?? . Taking the mean value of the path-wise strong error yields \[\epsilon^{\text{Sim}}(\Delta t_R)=1/U\sum_{u=1}^{U}e^{\text{Sim}}(u,\Delta t_R).\]

The comparison of the convergence order of fEMM and method ?? is shown in 1.

2 shows the density (PDF) 43 of the spectral distribution of the exact solution \(X_t\) of 42 recovered from it’s Cauchy transform at two different time points and their approximation calculated by ?? in \(\mathcal{M}^{500}(\mathbb{R})\). The bars show an estimation of the PDF via the eigenvalues of the numerical solution \(\overline{X}_t\), calculated by ?? at different time points. The line is the PDF of the exact solution \(X_t\).

8.2 Geometric Brownian Motion I↩︎

Let’s consider the equation \[\label{ex:geo1} dX_t= \theta X_tdt + \sqrt{X_t}dW_t\sqrt{X_t}\tag{44}\] with \(X_0=I\). This examples demonstrates numerically the order \(\gamma=1\) of the method developed in 23. In every single time step, the spectral properties of \(\overline{X}_k\) has to be computed in order to calculate the single and double operator integrals 28 . The difficulty of this equations is, that if \(\theta\) is to small, the eigenvalues turn complex. To confirm \(\gamma=1\) we choose \(\theta\) large enough. The equations is solved in the von Neumann algebra \(\mathcal{M}^5(\mathbb{R})_{sa}\). 3 shows the order of convergence.

The reference solution \(\overline{X}_T(u,\Delta t_0)\) was calculated by \(\Delta t_0=1/2^{12}\). Further calculations were performed with \(\Delta t_R=1/2^{R+7}\) for \(R=1,2,3\). We then calculated the error \(e^{\text{Sim}}(u,\Delta t_R)=\|\overline{X}_T(u,\Delta t_R) - \overline{X}_T(u,\Delta t_0)\|\) for each path number \(u\). As in the previous example, we calculate the mean value of the pathwise error \(e^{\text{Sim}}(u,\Delta t_R)\). The plot in 3 shows the expected order of convergence.

8.3 Geometric Brownian Motion II↩︎

Consider the case \(d>1\) with linear diffusion coefficients of the form \[\label{ex:geo2} dX_t= \theta X_tdt + X_tdW_t + dW_t X_t\tag{45}\] with start value \(X_0=I\). For analytical insights to the spectral distribution of the solution \(X_t\) we again refer to kargin?. We apply method fSM ?? and set \(N=10\), \(\theta=1\). The subdivision method fSM is realized by taking time step values \(\Delta t_R = 2^R\cdot 10^{-3}\) for \(R=0,1,2,3\) and \(\delta t_R=\Delta t_R^2\). We adopt the abbreviations from the first example. The underlying path of the Brownian motion was calculated by a time step \(10^{-6}\) (\(R=0\)). Out of this path, the necessary increments of the Brownian motion for the application of fSM are calculated in dependence of the time step \(\Delta t_R\). As in the previous example, we consider \(U\in\mathbb{N}\) different paths. Due to the computational complexity of fSM we perform an a posterior error calculation. For each path \(u\in\mathbb{N}\) we calculate \[\begin{gather} e^{GeoII}(u,\Delta t_R)=\|\overline{X}_T(u,\Delta t_R)-\overline{X}_T(u,\Delta t_{R-1})\|_1 =\\= \frac{1}{N}\mathbb{E}\left(\text{tr}(|\overline{X}_T(u,\Delta t_R)-\overline{X}_T(u,\Delta t_{R-1})|)\right). \end{gather}\]\(\overline{X}_T(u,\Delta t_R)\) denotes the numerical solution by an appropriate method with path \(u\) and time step \(\Delta t_R\) at time point \(T\). We consider \(U=2000\) number of different paths and approximate the posterior error as the mean value \(\epsilon^{\text{GeoII}}(\Delta t_R)=\frac{1}{U}\sum_{u=1}^U e^{\text{GeoII}}(u,\Delta t_R)\). The estimation of the order of convergence is then calculated by \[\begin{align} \label{eq:aposterioriestim} \gamma \approx \overline{\gamma}(\Delta t_R)=\frac{\log\left(\frac{\epsilon^{\text{GeoII}}(\Delta t_R)}{\epsilon^{\text{GeoII}}(\Delta t_{R-1})}\right)}{\log(2)}. \end{align}\tag{46}\] The results are listed in 1 and visualized in 4. The simulation shows that \(\overline{\gamma}(\Delta t_3)\approx 0.961\) by fSM and \(\overline{\gamma}(\Delta t_3)\approx 0.501\) by fEMM. The results confirm the expected strong convergence orders of \(\gamma=0.5\) for fEMM and \(\gamma=1\) for fSM.

 
 

8.4 CIR Equation↩︎

Consider the following nonlinear case with \(d>1\), \[\label{eq:CIR} dX_t= (a-bX_t)dt + \frac{\sigma}{2}\sqrt{X_t}dW_t + dW_t\frac{\sigma}{2}\sqrt{X_t}\tag{47}\] with start value \(X_0=I\). For analytical insights to the spectral distribution and the existence of the solution \(X_t\) we refer to [14]. We perform the simulation with \(N=2\), \(a=I\), \(b=0.1\), \(\sigma=0.2\). Similar to the example in 8.3, we proceed by applying method fSM ?? with time steps \(\Delta t_R = 2^R\cdot 10^{-3}\) for \(R=0,1,2,3\) and \(\delta t_R=\Delta t_R^2\). The underlying path of the Brownian motion was calculated by a time step \(10^{-6}\). The error for path \(u\in\mathbb{N}\) is calculated at \[\begin{gather} e^{\text{CIR}}(u,\Delta t_R)=\|\overline{X}_T(u,\Delta t_R)-\overline{X}_T(u,\Delta t_{R-1})\|_1 =\\= \frac{1}{N}\mathbb{E}\left(\text{tr}(|\overline{X}_T(u,\Delta t_R)-\overline{X}_T(u,\Delta t_{R-1})|)\right). \end{gather}\] \(\overline{X}_T(u,\Delta t_R)\) denotes the numerical solution by fEMM or fSM with path \(u\) and time step \(\Delta t_R\) at time point \(T\). We consider \(U=2000\) number of different paths and approximate the a posteriori error as \(\epsilon^{\text{CIR}}(\Delta t_R)=\frac{1}{U}\sum_{u=1}^Ue^{\text{CIR}}(u,\Delta t_R)\). The convergence order is numerically estimated by \[\begin{align} \label{eq:aposterioriestimCIR} \gamma \approx \overline{\gamma}(\Delta t_R)=\frac{\log\left(\frac{\epsilon ^{\text{CIR}}(\Delta t_R)}{\epsilon^{\text{CIR}}(\Delta t_{R-1})}\right)}{\log(2)}. \end{align}\tag{48}\] The results are listed in 2 and visualized in 5. The simulation show that \(\overline{\gamma}(\Delta t_3)\approx 0.999\) by fSM and \(\overline{\gamma}(\Delta t_3)\approx 0.523\) by fEMM.

 
 

References↩︎

---

1. Correspondig author email: michael.wibmer@hm.edu↩︎