I have a question about the order of convergence of the Euler-Maruyama scheme and how one verifes this numerically.
I have read that the Euler-Maruyama scheme is mean-squared convergent of order 1/2 where (I think) mean-square convergence is defined by:
$$ \max_{n=1\dots N} \lVert X(t_n) -Y_n \rVert_{L_2} \leq C h^{1/2}, $$ where $ \lVert Z \rVert_{L_2} = \left( \mathbb{E}|Z|^2\right)^{1/2}$, $X(t_n)$ is the solution to an SDE at timepoint $t_n$ and $Y_n$ is the Euler-Maruyama approximation at $t_n$.
I want to verify this numerically and my approach is to approximate the norm with a Monte Carlo simulation by simulating $M$ sample paths and computing
$$ \lVert X(t_n) -Y_n \rVert_{L_2} \approx \left( \frac{1}{M} \sum_{j=1}^M |X^j(t_n) -Y^j_n|^2 \right)^{1/2} $$
($j$ is the number of the simulated sample path) for a sequence of different grids and for each grid take the maximum of the Monte Carlo means. Instead of using the true solution to the SDE I calculate the Euler-Maruyama solution on a very fine grid and use this as the reference solution. I then simulate the EM-approximation on coarser grids (but using the driver process that I used for the reference solution). This however gives me strange results (the convergence plot doesn't indicate any convergence).
It seems to work (I get order of convergence $1/2$) if I instead take the $|X^j(T) - Y^j_N|$ (the absolute value of the difference at the end point for each sampel path $j$) and then take the avergage of these for each grid. My question is if I have misunderstood something or my approach is incorrect?
Thank you and kind regards,