Discretization formula for system of differential equations. “Solution to one of these is the initial condition of the other”. In which sense?

Consider the following stochastic differential equation $$$$dy=\left(A-\left(A+B\right)y\right)dt+C\sqrt{y\left(1-y\right)}dW\tag{1}$$$$ where $$A$$, $$B$$ and $$C$$ are parameters and $$dW$$ is a Wiener increment.
Equation $$(1)$$ will be our point of reference in what follows.

Now, first let us consider a "method" for equation $$\left(1 \right)$$ which can be described by the following one-step discretization scheme: $$$$y_{n+1}=y_n+\left(A-\left(A+B\right)y_n\right)\Delta t +C\sqrt{y_n\left(1-y_n\right)}\Delta W_n + D\left(y_n\right)\left(y_n-y_{n+1}\right)\tag{2}$$$$ where $$\Delta t$$ is the length of the time discretization interval, $$\Delta W_n$$ is a Wiener increment and $$D(y_n)$$ is the system of control functions and takes the form $$D(y_n)=d^0(y_n)\Delta t + d^1\left(y_n\right)|\Delta W_n|$$ with $$d^1(y)= \begin{cases} C\sqrt{\frac{1-\varepsilon}{\varepsilon}}\hspace{0.5cm}\text{if }y<\varepsilon\\ C\sqrt{\frac{1-y}{y}}\hspace{0.5cm}\text{if }\varepsilon\le y<\frac{1}{2}\\ C\sqrt{\frac{y}{1-y}}\hspace{0.5cm}\text{if }\frac{1}{2}\le y\le 1-\varepsilon\\ C\sqrt{\frac{1-\varepsilon}{\varepsilon}}\hspace{0.5cm}\text{if }y>1-\varepsilon \end{cases}$$ At this point, let us consider a "method" which decomposes $$\left(1\right)$$ into two equations. Specifically, the first equation is a stochastic one, that consists of the diffusion term of $$\left(1\right)$$ only (see eqtn $$\left(3\right)$$), while the second one is an ordinary differential equation (see eqtn $$\left(4\right)$$) that consists of the drift part of $$\left(1\right)$$. We have:

$$$$dy_1=C\sqrt{y_1\left(1-y_1\right)}dW\tag{3}$$$$ $$$$dy_2=\left(A-\left(A+B\right)y_2\right)dt\tag{4}$$$$

This last method approximates the solution to $$\left(3\right)$$ at each time step using $$\left(2\right)$$ (and numerical solution to $$\left(3\right)$$ is used as the initial condition in $$\left(4\right)$$), while $$\left(4\right)$$ can be solved using the Euler method. Thus, such a method can be described by the following one step discretization formula: $$y_{n+1}=y_n+\left(A-\left(A+B\right)y_n\right)\Delta t + \dfrac{C\sqrt{y_n\left(1-y_n\right)}\Delta W_n}{1+d^1\left(y_n\right)|\Delta W_n|}\left(1-\left(A+B\right)\Delta t\right)\tag{5}$$

My doubts:

1. I cannot understand in which way the last method approximates solution to $$\left(3\right)$$ at each time step using $$\left(2\right)$$. Could you please explicit such an approximation? How is it obtained by means of $$\left(2\right)$$?
2. In which sense numerical solution to $$\left(3\right)$$ is used as the initial condition in $$\left(4\right)$$? Which is such an initial condition?
3. Could you please explicit the way in which solution to $$\left(3\right)$$ and solution to $$\left(4\right)$$ are combined so as to obtain discretization formula $$\left(5\right)$$?
• Does your source have any remark on how this method improves Euler-Maruyama, Milshtein and possibly its derivative-free variant? Is the construction of $d^1y$ related to the singular derivative of this variation coefficient? – Lutz Lehmann Jul 23 '20 at 11:31
• In summary, in my source it is stated that methods such as Euler-Maruyama and Milstein are unable to preserve the boundaries to the solution of Wright-Fisher SDE (which must be $[0,1]$), while the above method shows to be able to preserve such boundaries. It is not explicitly stated whether construction of $d^1y$ is related to the singular derivative of this variation coefficient. Thank you a lot if you could help me, it would be very very appreciated @LutzLehmann – Strictly_increasing Jul 23 '20 at 11:57
• That appears as a reasonable goal, as $\Delta W_t=W_{t+\Delta t}-W_t$ is unbounded over the full probability space. Are any claims made on the order of the resulting method, I think preserving the order 0.5 of Euler-Maruyama can be obtained with easier modifications, simple cut-offs? – Lutz Lehmann Jul 23 '20 at 12:03
• $ΔW_t(\omega)$ as random variable for fixed $t$ and $Δt$ is unbounded (also in the essential sense) as a normally distributed variable. The pure Euler update is thus likewise unbounded, leaving the interval $[0,1]$ with some positive probability. The fraction in the final formula makes this contribution bounded, but it is not immediately visible how that results in values inside the interval. – Lutz Lehmann Jul 23 '20 at 16:02
• This topic might also get contributions on scientific computations scicomp.SE or on quantitative fincance, quant.SE. – Lutz Lehmann Jul 23 '20 at 18:53