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.
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.