Consider the following stochastic differential equation \begin{equation} dy=\left(A-\left(A+B\right)y\right)dt+C\sqrt{y\left(1-y\right)}dW\tag{1} \end{equation} 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: \begin{equation} 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} \end{equation} 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:

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

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)$?
  • 1
    $\begingroup$ 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? $\endgroup$ Jul 23, 2020 at 11:31
  • $\begingroup$ 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 $\endgroup$ Jul 23, 2020 at 11:57
  • 1
    $\begingroup$ 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? $\endgroup$ Jul 23, 2020 at 12:03
  • 1
    $\begingroup$ $Δ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. $\endgroup$ Jul 23, 2020 at 16:02
  • 1
    $\begingroup$ This topic might also get contributions on scientific computations scicomp.SE or on quantitative fincance, quant.SE. $\endgroup$ Jul 23, 2020 at 18:53


You must log in to answer this question.