# 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? Jul 23, 2020 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 Jul 23, 2020 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? Jul 23, 2020 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. Jul 23, 2020 at 16:02
• This topic might also get contributions on scientific computations scicomp.SE or on quantitative fincance, quant.SE. Jul 23, 2020 at 18:53