Evaluating multiple integrals after Itô-Taylor expansion Consider an autonomous, scalar stochastic differential equation (SDE):
$$ d[X(t)] = f[X(t)]\textrm{d}t + g[X(t)]\textrm{d}W(t) $$
Consider also a scalar function $U[X(t)]$ of the solution of the SDE.
The integral form of the Itô formula is:
$$ U[X(t)] = U[X(t_0)] + \int_{t_0}^{t} L^0 U[X(s)] \textrm{d}s + \int_{t_0}^{t} L^1 U[X(s)] \textrm{d} W(s) $$
With the differential operators:
$$ L^0 = f \frac{\partial}{\partial X} + \frac{1}{2} g^2 \frac{\partial^2}{\partial X^2} $$
$$ L^1 = g \frac{\partial}{\partial X} $$
We consider a third-order expansion for the case of a SDE with a constant function $g[X(t)]$. After several applications of the Itô formula to the integrands $U[X(t)] = f[X(t)]$ and $U[X(t)] = g[X(t)]$ in the SDE, we obtain:
$$ \begin{align}
X\left(t\right) = X\left(t_0\right) & + fI_0 + gI_1 + \left(ff^\prime+\frac{1}{2}g^2f^{\prime\prime}\right)I_{00} + \left(gf^\prime\right)I_{10} \\ 
& +\left[f\left(ff^{\prime\prime}+f^\prime f^\prime+\frac{1}{2}g^2f^{\prime\prime\prime}\right)+\frac{1}{2}g^2\left(ff^{\prime\prime\prime}+3f^\prime f^{\prime\prime}+\frac{1}{2}g^2f^{\prime\prime\prime\prime}\right)\right]I_{000} \\ 
& +\left[g\left(ff^{\prime\prime}+f^\prime f^\prime+\frac{1}{2}g^2f^{\prime\prime\prime}\right)\right]I_{100} \\ 
& +\left[f\left(gf^{\prime\prime}\right)+\frac{1}{2}g^2\left({gf}^{\prime\prime\prime}\right)\right]I_{010} \\ 
& +\left[g\left(gf^{\prime\prime}\right)\right]I_{110} \\ 
& + R 
\end{align} $$
Where $f$, $g$, and the derivatives of $f$ are evaluated at $X(t_0)$.
This equation includes several stochastic and deterministic integrals. Some are straightforward to evaluate:
$$ \begin{align}
I_0 & = \int_{t_0}^{t}{\mathrm{d}s} = \Delta t \\
I_1 & = \int_{t_0}^{t}{\mathrm{d}W(s)} = \Delta W\left(t\right) \\
I_{00} & = \int_{t_0}^{t}\int_{t_0}^{s}{\mathrm{d}u\mathrm{d}s} = \frac{1}{2}{(t-t_0)}^2 = \frac{{\Delta t}^2}{2} \\
I_{000} & = \int_{t_0}^{t}\int_{t_0}^{s}\int_{t_0}^{u}{\mathrm{d}v\mathrm{d}u\mathrm{d}s}=\frac{{\Delta t}^3}{6}
\end{align} $$
But others, not so much. They are the double or triple integrals of an integrand equal to 1 with respect to either the time $\mathrm{d} s$ or a Brownian motion $\mathrm{d} W(s)$:
$$ \begin{align}
I_{10} & = \int_{t_0}^{t}{\int_{t_0}^{s}{\mathrm{d}W(u)}\mathrm{d} s} \\
I_{100} & =\int_{t_0}^{t}\int_{t_0}^{s}\int_{t_0}^{u}{\mathrm{d}W(v)\mathrm{d}u\mathrm{d}s} \\
I_{010} & =\int_{t_0}^{t}\int_{t_0}^{s}\int_{t_0}^{u}{\mathrm{d}v\mathrm{d}W(u)\mathrm{d}s} \\
I_{110} & =\int_{t_0}^{t}\int_{t_0}^{s}\int_{t_0}^{u}{\mathrm{d}W(v)\mathrm{d}W(u)\mathrm{d}s}
\end{align} $$
Could somenone please explain how to evaluate these last four integrals?
My hope is to obtain a scalar, for instance $\sqrt{\Delta t}N(0,1)$ as in $I_1$, in order to use Itô-Taylor expansion as a numerical integration scheme.
Many thanks
 A: You can proceed in the same manner as for the classical nested integrals. For instance, the first one is given by
$$
I_{10} = \int_{t_0}^t\int_{t_0}^s\mathrm{d}W_u\mathrm{d}s = \int_{t_0}^t(W_s-W_{t_0})\mathrm{d}s = \int_{t_0}^tW_s\mathrm{d}s - (t-t_0)W_{t_0}
$$
Thanks to Itô's lemma applied to the quantity $tW_t$, it is possible to rewrite the remaining integral as follows :
$$
\int_{t_0}^tW_s\mathrm{d}s = \int_{t_0}^t\mathrm{d}(sW_s) - \int_{t_0}^ts\mathrm{d}W_s = tW_t-t_0W_{t_0} - \int_{t_0}^ts\mathrm{d}W_s
$$
hence
$$
I_{10} = t(W_t-W_{t_0}) - \int_{t_0}^ts\mathrm{d}W_s
$$
Recalling that $\int_{t_0}^ts\mathrm{d}W_s$ is gaussian with
$$
\mathrm{Var}\left[\int_{t_0}^ts\mathrm{d}W_s\right] = \int_{t_0}^ts^2\mathrm{d}s = \frac{1}{3}(t^3-t_0^3)
$$
by Itô's isometry, we can conclude that
$$
I_{10} \sim \mathcal{N}\left(0,\frac{2}{3}t^3+t^2t_0+\frac{1}{3}t_0^3\right)
$$
since $\lambda_1\mathcal{N}(\mu_1,\sigma_1^2) + \lambda_2\mathcal{N}(\mu_2,\sigma_2^2) = \mathcal{N}(\lambda_1\mu_1+\lambda_2\mu_2,\lambda_1^2\sigma_1^2+\lambda_2^2\sigma_2^2)$ for independent gaussians. The other integrals can be handled in the same way.
Finally, note that, apart from the expansion you are trying to implement, numerical methods based on finite elements are also possible for stochastic differential equations (see e.g. here or here).

Addendum.
Here are further hints for the main obstacles with the other integrals. Note that we have :
$$
I_{100} = \int_{t_0}^tI_{10}(s)\mathrm{d}s = \int_{t_0}^t(sW_s-W_{t_0})\mathrm{d}s - \int_{t_0}^t\int_{t_0}^s\sigma\mathrm{d}W_\sigma\mathrm{d}s
$$
The integral containing the term $sW_s$ can be tackled by the same procedure as for $I_{10}$, that is with the help of Itô's lemma for the quantity $s^2W_s$. But what to do with $X_t := \int_{t_0}^t\int_{t_0}^s\sigma\mathrm{d}W_\sigma\mathrm{d}s$ (this kind of integral will also appear in $I_{010}$) ? Recall that $\int_{t_0}^s\sigma\mathrm{d}W_\sigma = W_{\phi(s)}$, with $\phi(s) := \frac{1}{3}(s^3-t_0^3)$, hence thanks to the change of time $\tau = \phi(s)$ :
$$
X_t = \int_{t_0}^tW_{\phi(s)}\mathrm{d}s = \int_{\phi(t_0)}^{\phi(t)}W_\tau\frac{\mathrm{d}\tau}{\dot{\phi}(\phi^{-1}(\tau))} = \int_0^{\phi(t)}\frac{W_\tau}{(3\tau+t_0^3)^{3/2}}\mathrm{d}\tau
$$
Applying Itô's lemma to $\psi(\tau)W_\tau$, where $\psi(\tau) := \frac{-2/3}{\sqrt{3\tau+t_0^3}}$, leads to
$$
X_t = \int_0^{\phi(t)}\mathrm{d}(\psi(\tau)W_\tau) - \int_0^{\phi(t)}\psi(\tau)\mathrm{d}W_\tau \sim \mathcal{N}\left(0,\phi(t)\psi(\phi(t))^2 + \int_0^{\phi(t)}\psi(\tau)^2\mathrm{d}\tau\right),
$$
which I let you compute explicitly.
As for the last integral, namely $I_{110}$, it contains the following quantity :
$$
\int_{t_0}^t\mathrm{d}W_s\int_{t_0}^s\mathrm{d}W_\sigma = \int_{t_0}^t(W_s-W_{t_0})\mathrm{d}W_s = \frac{1}{2}(W_t^2-W_{t_0}^2)-\frac{1}{2}(t-t_0) + W_{t_0}(W_t-W_{t_0}),
$$
where $\int_{t_0}^tW_s\mathrm{d}W_s = \frac{1}{2}(W_t^2-W_{t_0}^2)-\frac{1}{2}(t-t_0)$ is a well-known result obtained by Itô's lemma on $W_s^2$; unfortunately, it is not gaussian. The last term is not gaussian either (unless $t_0 = 0$, because $W_0 = 0$); it is to be noticed that
$$
\mathbb{E}[W_{t_0}(W_t-W_{t_0})] = \mathbb{E}\left[\int_{t_0}^tW_{t_0}\mathrm{d}W_s\right] = \int_{t_0}^t\mathbb{E}[W_{t_0}]\mathbb{E}[\mathrm{d}W_s] = 0
$$
and
$$
\mathrm{Var}[W_{t_0}(W_t-W_{t_0})] = \mathrm{Var}\left[\int_{t_0}^tW_{t_0}\mathrm{d}W_s\right] = \mathbb{E}\left[\left(\int_{t_0}^tW_{t_0}\mathrm{d}W_s\right)^2\right] = \int_{t_0}^tt_0\mathrm{d}s = t_0(t-t_0)
$$
