Stochastic integral formula numerical approximation dependence on $\Delta t$

Question

I have the equation $$\frac1n B_t^n=\int_0^t B_s^{n-1} dB_s+\frac{n-1}{2}\int_0^t B_s^{n-2} ds$$

for $n\geq 2$. Of course it works with $n=1$ also if we drop the last integral. I derived from applying the Ito formula to $d(\frac1n B_t^n)$ where $B_t$ is standard Brownian motion starting at zero.

I'm trying to approximate both sides of the equation numerically, and it generally works pretty well. The difference can sometimes be large when $n$ is large or $\Delta t$ isn't too small.

I'd like to know how the difference in the left and right sides of the equation depends on the step size $\Delta t$.

So I want to understand the difference $$D\equiv \sum_{i=1}^{N-1}B_{t_i}^{n-1}(B_{t_{i+1}}-B_{t_i})+\frac{n-1}2\sum_{i=1}^{N-1}B_{t_i}^{n-2}\Delta t - \frac1nB_{t_{N}}^n$$ where the $\{t_i\}_{i=0}^N$ partition the interval $(0,t)$ and $\Delta t =t_{i+1}-t_i.$ I start the summation at $i=1$ since $B_{t_0}=B_0=0.$

I thought that maybe we could understand it in terms of the approximation of $B_{t_{i+1}}-B_{t_i}$ as $\mathcal N(0,\Delta t)$ random variables, which is how I am numerically simulating this. I.e. let $X_{i}=B_{t_{i}}-B_{t_{i-1}}$ be iid $\mathcal N(0,\Delta t)$ for $i\in\{1,2,\cdots,N\}$, and we get

$$D=\sum_{i=1}^{N-1}X_{i+1}\left(\sum_{k=1}^{i}X_k\right)^{n-1} +\frac{n-1}2\Delta t\sum_{i=1}^{N-1}\left(\sum_{k=1}^{i}X_k\right)^{n-2} - \frac1n\left(\sum_{k=1}^{N}X_k\right)^n.$$

I took the expectation of this and used the fact that the $X_i$ are mean zero and iid, that $S_i=\sum_{k=1}^i X_k \sim \mathcal N(0,i\Delta t)$ and $t=N\Delta t$ to get $$\begin{aligned} E(D)&= \sum_{i=1}^{N-1}E(X_{i+1})E\left[S_i^{n-1}\right] +\frac{n-1}2\Delta t\sum_{i=1}^{N-1}E\left[S_i^{n-2}\right] - \frac1nE\left[S_N^n\right]\\ &=0+\frac{n-1}2\Delta t\sum_{i=1}^{N-1}m_{n-2}^{i\Delta t} - \frac1nm_{n}^{t}\\ \end{aligned}$$ where $m_n^{\sigma^2}$ is the $n^\text{th}$ moment of $\mathcal N(0,\sigma^2)$ distribution.

Using $m_n^{\sigma^2}=\sigma^n(n-1)!!$, this simplifies to

$$E(D)=\frac{n-1}2\Delta t\sum_{i=1}^{N-1}(n-3)!!{(i\Delta t)}^{(n-2)/2} - \frac1n(n-1)!!t^{n/2}$$ and then using $\Delta t=\frac t N$ again, it becomes $$E(D)=(n-1)!! \ t^{n/2}\left(\frac12N^{-n/2}\sum_{i=1}^{N-1}{i}^{n/2-1} - \frac1n\right).$$

Just to get a particular example of $E(D)$, I'll let $t=1$, $\Delta t=0.01$, $n=3$, and thus $N=100$ to get $$ \begin{aligned} E(D)&=0.001\sum_{i=1}^{99}{i}^{1/2} - \frac23 \cdot\\ &\approx 0.001\cdot\frac23\cdot 99\cdot\sqrt{99+\frac32} - \frac23\\ &\approx -0.005019. \end{aligned}$$ Note that I used an approximation for $\sum_{i=1}^N\sqrt{i}\approx \frac23N\sqrt{N+\frac32}$ that I found elsewhere here on MathSE. $E(D)$ appears to be close to $-\frac12\Delta t$ for $t=1$, $n=3$. The variance of $D$ (although I haven't calculated it) seems to be sufficiently large relative to its mean that I cannot quickly assess if $E(D)\approx-0.005019$ is matched by numerics. (EDIT: For these parameter values $\text{Var}(D)\approx0.01$ from numerics. However, my confidence interval estimates of $E(D)$ from numerics $E(D_{sim})\pm3\sqrt{\frac{\text{Var}(D_{sim})}{n}}$ almost never capture $-0.005019$ for $N=10,000$, so I may have an error somewhere, either in my numerical code or something analytical here.)

I'm not sure if calculating $E(D)$ is even useful or where to go from here. Also, maybe I have a big mistake or conceptual misunderstanding somewhere.

What I am hoping is that I can get an estimate of this difference in terms of how it depends on $\Delta t$ and maybe the particular $B_t$ realization. E.g. maybe $D$ is approximately normal with mean $\mu(\Delta t,t,n)$ and variance $\sigma^2(\Delta t,t,n)$ are some tractable functions of $\Delta t,t,n.$

---

EDIT: I fixed the approximation as per @ClaudeLeibovici answer-comment below. Simply by trial an error I found that $$\sum_{i=1}^{N-1}i^{\frac n2-1}\approx \frac2nN^\frac n2-\frac12N^{\frac n2-1}.$$ Plugging this into $E(D)$ gives $$E(D)\approx -\frac14(n-1)!! \ t^{\frac n2-1}\Delta t.$$ I have yet to attempt to calculate $\text{Var}(D)$ because I think it will probably be tedious. I did calculate $E(\frac1nB_t^n)=\frac1nt^\frac n2 (n-1)!!$ and $\text{Var}(\frac1nB_t^n)=\frac{1}{n^2}t^n(2n-1)!!$ but not sure if that is helpful.

Here is the code I'm using in Octave. Maybe there is an error in it.

n=3;
t=1;
dt=0.01; 
N=ceil(t/dt); 
X=randn(1,N)*sqrt(dt);

%approx integral of B^(n-1)dB
integral1=sum(cumsum(X(1:N-1)).^(n-1).*X(2:N))

%approx integral of B^(n-2)dt
integral2=sum(cumsum(X(1:N-1)).^(n-2))*dt;

D=integral1+(n-1)/2*integral2-1/n*sum(X(1:N))^n

Answer 1

This is not an answer but it is too long for a comment.

I do not know where on MSE you did find the approximation $$S_n=\sum_{i=1}^n\sqrt{i}\approx \frac23n\sqrt{n+\frac54}$$ but it seems that $$T_n=\sum_{i=1}^n\sqrt{i}\approx\frac23n\sqrt{n+\frac32}$$ would be a better approximation.

For example, for $n=1000$, the first formula would give $21095.023$ while the second gives $21097.657$, the "exact" value being $21097.456$.

If you are concerned with very large values of $n$, taking into account that $$U_n=\sum_{i=1}^n\sqrt{i}=H_n^{\left(-\frac{1}{2}\right)}$$ where appear the generalized harmonic numbers. Using the corresponding asymptotics, you would get $$U_n=\frac{2}{3} \sqrt{n} \left(n+\frac{3}{4}\right)+\zeta \left(-\frac{1}{2}\right)+O\left(\frac 1 { \sqrt n}\right)$$ which, for $n=1000$ gives $21097.455$.