Euler-Maruyama method is supposed to be an extension of the Euler method for ODE, but applied to SDE. This means that if we have an equation:
$$ dY_t = Y_t dW_t $$
where $W_t$ is the Wiener process, discretizing the total time T into $N$ intervals ($T = N \cdot h$), we could approximate the change in $Y$ from $t=0$ to $t=1$:
$$ d Y_1 = Y_1 - Y_0 \approx Y_1 \cdot (W_1 - W_0) $$
Now to evaluate this (implicit) equation for $Y_1$ I would have to generate a random number for $W_1 - W_0$ which is a normal distribution $~ \mathcal{N}(0,1)$, as in the definition of the Wiener process.
However in the wikipedia page (https://en.wikipedia.org/wiki/Euler%E2%80%93Maruyama_method) as well as in other sources the E-M method seems to dictate a random number from $\mathcal{N}(0,1)$ TIMES $\sqrt{h}$.
Why is the square root of the timestep $h$ there? I suspect it has something trivial to do with the standard deviation, but I wasn't able to derive this result, since the expactation value of $W_t - W_s$ is $0$. Would you help a fellow out?
First, the Ito interpretation prevents you from using backwards Euler that easily, you can only use forward Euler with such a literal translation.
Second, the difference you are considering is actually $W_{t_1}-W_{t_0}=W_h-W_0$, and the increments of the Wiener process follow the law $$ W_s-W_t\sim N(0,|s-t|)\sim N(0,1)·\sqrt{|s-t|} $$ which should explain where the square root comes from.
If one wanted to apply the implicit Euler method, note that while $(1+ΔW)$ is the factor from $Y_0$ to $Y_1$, its inverse in the Ito formalism computes as $(1-ΔW+ΔW^2+...)\simeq(1-ΔW+Δt)$ so that the resulting implicit Euler step reads as $$ Y_0=(1-ΔW+Δt)Y_1\iff Y_1=Y_0-Y_1Δt+Y_1ΔW $$ One needs to apply similar corrections for more complex SDE.