11

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)$?
Strictly_increasing
  • 555
  • 1
    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
  • 1
    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
  • It is stated that numerical tests suggest that the algorithm still converges strongly with order 1/2. When saying that $\Delta W_t$ is unbounded over the full probability space, you mean that it is not almost sure that it stands between $0$ and $1$, don't you? @LutzLehmann Do you have any suggestion/answer to my doubts above? – Strictly_increasing Jul 23 '20 at 12:48
  • 1
    $Δ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
  • Of course, and your final doubt is the same I have. I was curiuous about how to get exactly there (to that correction at denominator). Anyway, thank you a lot for all your help so far, something is gonna be clearer now. If you have some suggestion on how to proceed so as to try to understand the mathematical derivation of that denominator correction, I would be extremely glad :) @LutzLehmann – Strictly_increasing Jul 23 '20 at 16:15
  • Pardon do you think that I could receive more feedbacks if I posted such a question in some other Stack Exchange environment? Perhaps some more "numerical-methods-oriented" enviroment? Does it exist at all? @LutzLehmann – Strictly_increasing Jul 23 '20 at 17:19
  • 1
    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

0 Answers0