Show that the Runge-Kutta method for 1-d SDE is weakly consistent with order 2 
Consider the following Runge-Kutta methods for one-dimensional SDE $dX_t=b(x,t)dt+\sigma(x,t)dW_t$:
\begin{multline}
X_h^{k+1} = X_h^k+\frac{1}{2}(b(X_h^k)+b(Y))h
  +\frac{1}{4}\left[\sigma(Y^+)+\sigma(Y^-)+2\sigma(X_h^k)\right]\xi^{k+1}
 \\ +\frac{1}{4\sqrt{h}}\left[\sigma(Y^+)-\sigma(Y^-)\right]\left[(\xi^{k+1})^2-h\right]\tag1
\end{multline}
where $$Y=X_h^k+hb(X_h^k)+\sigma(X_h^k)\xi^{k+1}\tag{2a}$$ and $$Y^\pm=X_h^k+hb(X_h^k)\pm\sigma(X_h^k)\sqrt{h}\tag{2b}$$ and $\xi$ are noises with $\mathbb{E}[\xi]=0$ and $\mathbb{E}[\xi^2]=h$.

I'd like to show that when $b$ is smooth (with bounded derivatives) and $\sigma$ is constant, the Runge-Kutta method is weakly consistent with order $p=2$, i.e. it satisfies the following bound with $p=2$,
$$
\left|\mathbb{E}\left[(F_h(k,x,\xi^{k+1})^m\right]
 -\mathbb{E}\left[(\Delta_{kh}^{(k+1)h}X)^m|X^{(kh)} = x\right]\right|
\leq Kh^{p+1},~~~\forall m\leq 2p+1  \tag3
$$
My way to justify this is to take the first and second moment of $F(x,\xi)$ and $\Delta_0^hX$, where
\begin{multline}
F(x,\xi) = \frac{1}{2}(b(X_n^k)+b(Y))h
  +\frac{1}{4}\left[\sigma(Y^+)+\sigma(Y^-)+2\sigma(X_h^k)\right]\xi^{k+1}
 \\ +\frac{1}{4\sqrt{h}}\left[\sigma(Y^+)-\sigma(Y^-)\right]\left[(\xi^{k+1})^2-h\right]  \tag{4a}
\end{multline}
and
\begin{multline}
\Delta_0^hX = b(x)h + \int_0^h\int_0^s\mathcal{L}b(X^r)drds 
  + \int_0^h\int_0^sb'(X^r)\sigma(X^r)dW_rds + \sigma\Delta_0^h W
  \\ + \int_0^h\int_0^s\mathcal{L}\sigma(X^r)drdW_s+\int_0^h\int_0^s\sigma'(X^s)\sigma(X^r)dW_rdW_s  \tag{4b}
\end{multline}
Note that the $drds$ term is $O(h^2)$, $dW_rds$ and $drdW_s$ terms are $O(h^{3/2})$, and $dW_rdW_s$ term is $O(h)$. Since $\sigma$ is constant, this reduces to
$$F(x,\xi)=\frac{1}{2}(b(X_n^k)+b(Y))h+\frac{1}{4}\left[\sigma(Y^+)+\sigma(Y^-)+2\sigma(X_h^k)\right]\xi^{k+1}  \tag{5a}$$
and
\begin{multline}
\Delta_0^hX = b(x)h+\int_0^h\int_0^s\mathcal{L}b(X^r)drds
  + \int_0^h\int_0^sb'(X^r)\sigma(X^r)dW_rds + \sigma\Delta_0^h W
  \\ + \int_0^h\int_0^s\mathcal{L}\sigma(X^r)drdW_s  \tag{5b}
\end{multline}
The first and second moments are
$$\mathbb{E}[F(x,\xi)]=\frac{b(X_h^k)+b(Y)}{2}h\tag{6a}$$
$$\mathbb{E}[\Delta_0^h X]=b(x)h+O(h^2)\tag{6b}$$
$$\mathbb{E}[F^2(x,\xi)]=\frac{h^2}{4}[b(X_h^k)+b(Y)]^2+\sigma^2\tag{6c}$$
$$\mathbb{E}[(\Delta_0^hX)^2]=h\sigma(x)+O(h^2)\tag{6d}$$
However, these does not satisfy the bound with $p=2$. Can anyone help me check where am I wrong? Thank you.
 A: Reduction of the method under the given assumption
Using the assumption $σ(x)=σ=const.$ simplifies the method to
\begin{align}
Y&=X_h^n+σ\Delta W_{nh}+b(X_h^n)h\\
X^{n+1}_h&=X_h^n+σ\Delta W_{nh}+\tfrac12[b(X_h^n)+b(Y)]h
\end{align}
Local truncation error and homogeneous components
Let's only analyze one step, the first step. If the exact solution $X_t$ is decomposed into components $X_t=x+X^{[1]}_t+X^{[2]}_t+...$ with $X^{[k]}_t\sim t^{k/2}$, then any (sufficiently smooth) function of $X_t$ can likewise be decomposed into homogeneous components,
\begin{align}
f(X_t)&=f(x)+F^{[1]}_t+F^{[2]}_t+...
\\
F^{[1]}_t&=f'(x)X^{[1]}_t\\
F^{[2]}_t&=f'(x)X^{[2]}_t+\tfrac12f''(x)(X^{[1]}_t)^2
\\
F^{[3]}_t&=f'(x)X^{[3]}_t+f''(x)X^{[1]}_tX^{[2]}_t+\tfrac16f'''(x)(X^{[1]}_t)^3
\\
F^{[4]}_t&=f'(x)X^{[4]}_t+f''(x)[2X^{[1]}_tX^{[3]}_t+(X^{[2]}_t)^2]
+\tfrac12f'''(x)(X^{[1]}_t)^2X^{[2]}_t+\tfrac1{24}f^{(4)}(x)(X^{[1]}_t)^4
\end{align}
For weak order $2$ the expectations of these terms $F^{[k]}$ have to match between the decompositions of exact and numerical solution for for $k=1,...,5$ so that the local errors $O(h^3)$ sum up to a global error $O(h^2)$. These terms have to match for any function $f$, so that the co-factors of the derivatives of $f$ have to match individually. The odd order terms $k=1,3,5$ are zero for symmetry reasons, this is true for the ODE as well as the symmetrically constructed method. If the first two components $X^{[k]}$, $k=1,2$ match exactly by construction of the method, then the only remaining expressions to match are the expectations of $X^{[1]}X^{[3]}$ and $X^{[4]}$.
Decomposition of the exact SDE solution
One can extract the homogeneous components from the integral form of the SDE, as time integration increases the half-step order by $2$
\begin{align}
X_t&=x+\int_0^tb(X_s)\,ds+σW_t\\
X^{[1]}_t&=σW_t\\
X^{[k]}_t&=\int_0^tB^{[k-2]}_s\,ds, ~~~ k=2,3,...
\\
\hline
X^{[2]}_t&=b(x)t\\
X^{[3]}_t&=\int_0^tb'(x)σW_s\,ds\\
X^{[4]}_t&=\int_0^t[b'(x)b(x)s+\tfrac12b''(x)σ^2W_s^2]\,ds
\end{align}
so that
\begin{align}
\Bbb{E}[X^{[1]}_tX^{[3]}_t]&=\int_0^tb'(x)σ^2s\,ds=\tfrac12b'(x)σ^2t^2
\\
\Bbb{E}[X^{[4]}_t]&=\int_0^t[b'(x)b(x)+\tfrac12b''(x)σ^2]s\,ds=\tfrac14[2b'(x)b(x)+b''(x)σ^2]t^2
\end{align}
Decomposition of the numerical step
Change notation to $\hat X_h=X^1_h$. Then $Y=x+\Delta Y$ and
\begin{align}
\hat X_t&=x+σW_t+b(x)t+\tfrac12[b'(x)ΔY+\tfrac12b''(x)(ΔY)^2+...]t
\\
&=x+σW_t+b(x)t+\tfrac12b'(x)σtW_t+\tfrac12[b'(x)b(x)t^2+\tfrac12b''(x)σ^2t(W_t)^2]
\end{align}
The required expectations match exactly.
