Let $ \frac{1}{2} < H < 1$ and let $B^H_t$ be a fractional Brownian motion with Hurst parameter $H$. Then the following stochastic differential equation $$\mathrm{d}Y_t = 5Y_t\mathrm{d}B^H_t, \quad Y_0 = 1$$ has as unique solution $Y_t = \exp(5B^H_t)$ in the sense of the Young (1936) integral.
Recently there has been some literature on solving such equations numerically. When considering the Euler-Maruyama scheme: \begin{align*} &Y_0^{(n)} = \xi \\ &Y_{(k+1)/n}^{(n)} = Y_{k/n}^n + 5Y_{k/n}^n\left( B^H_{(k+1)/n} - B^H_{k/n}\right), \quad k \in \{ 0,\dots, n-1\} \end{align*} on the interval $[0,1]$ and $n$ is the number of gridpoints. It has been proven (p. 5) that the following convergence result holds: \begin{equation*} n^{2H - 1} \| Y^{(n)} - Y \|_{\infty} \overset{\text{a.s.}}{\to} \dfrac{1}{2}\left(\sup_{t \in [0, 1]}\left\lvert\int_0^1 D_s Y_t \mathrm{d}s\right\rvert\right)\label{error1} \end{equation*} Where $D_sY_t$ is the Malliavin derivative of $Y_t$.
I want to empirically validate this result (and for a different, better approximation scheme). I have implemented the numerical scheme in R and got satisfactory results, but I am now stuck on how to exactly test convergence and how to compare my approximation with the exact result.
Is there any literature that details how to do this?
More specifically: (1) How can you compare the approximation with the exact result, given that it depends on the sample of the fractional Brownian motion and (2) How do you test the convergence result given above.
When dealing with normal ODEs, these questions are easy. But at the moment I am stuck...
Edit1: In the end I want to produce graphs similar to the ones given in this presentation on slides 14 through 17.