Testing Wiener process splitting in adaptive-step SDE integrators

Question

I am investigating various methods for adaptive-step integration of stochastic differential equations and trying to implement them. All of the papers that I've seen (e.g. H. Lamba, J. Comp. App. Math. 161(2), 417–430 (2003)) mention in some form that one has to always preserve the sampled values of the Wiener process. In other words, if some steps were rejected, one still has to go through all the sampled values that are located later in time.

For example, suppose that you are at time $t_0$ and tried to step to $t_1$, sampling the Wiener increment $dW_1$ for that. The step turned out to have an unacceptable error and was rejected. Now you try to make a step to $t_2 < t_1$. The proper procedure, it appears, is to sample an intermediate increment $dW_2$ for the step from $t_0$ to $t_2$ based on $dW_1$ (the exact formula is not important here), and take the increment for the step from $t_2$ to $t_1$ to be $dW_1 - dW_2$.

Now suppose that this step was rejected as well. Repeating the procedure above, you now have three increments, $dW_3$, $dW_2 - dW_3$ and $dW_1 - dW_2$ for the steps from $t_0$ to $t_3$, then to $t_2$ and then to $t_1$.

This part I have not seen explained, although it seems that all the authors imply that if your step to $t_3$ succeeds, your next step must not go past $t_2$ and should use the sampled increment $dW_2$ for that part. If you follow that logic, in the implementation of the integrator you will have to maintain a dynamic list of sampled increments, which increases the memory footprint (in an unpredictable way, too). This supposedly avoids bias in the sampled Wiener process. Let's call it approach A.

I tried to replace it with a simpler logic, where only the furthest sample in time ($dW_1$) is preserved. This way, in the situation above, after you sample $dW_3$, you forget $dW_2$, and for your next step the sampled increment only depends on $dW_1$. Let's call it approach B.

Sorry for the long explanation, now to my question. Supposedly my simplified logic should introduce bias in the solutions. The problem is that I haven't been able to detect it in any of the tests that I tried. Do you have any ideas of how to design a test that will be able to distinguish approach A from approach B?

Answer 1

I discuss the method you describe in more detail in this paper (Rackauckas and Nie 2017) as RSwM2. In that paper I am ever so slightly able to detect that it's sometimes doing something wrong, but since it only has issues with re-rejections it isn't that big of a deal. Those 3 methods (RSwM1, RSwM2, RSwM3) are now the basis of DiffEqNoiseProcess.jl and subsequently are the methods which build the adaptive timestepping in DifferentialEquations.jl.

Needless to say, I have a lot better tests now than I did back when the paper was published. The first thing to note is that this is a lot less of an issue if you switch to using PI-controlled adaptive timestepping (see Hairer II for details) as this reduces the chance of re-rejection tremendously. Also tweaking the adaptivity parameters helps as well. But, a very good test is to use adaptive timestepping to recapitulate the properties of the Brownian Bridge. Using the notation from the paper, if you start by putting a value into the future information stack $S_1$ at the final timepoint then you have to hit that value. So this means that starting with a value in $S_1$ should make the noise process be a Brownian Bridge. Since you know properties like the mean and variance of the Brownian Bridge at every point in time, it's easy to run a large number of simulations Monte Carlo and test for these properties in a very sensitive way. I have this test in my continuous integration testing suite and have found that it is super sensitive to getting the algorithm correct. Just to check, I ran it with the RSwM2 implementation:

W = BrownianBridge(0.0,1.0,0.0,1.0,0.0,0.0;rswm = RSWM(adaptivealg=:RSwM2))
prob = NoiseProblem(W,(0.0,1.0))
monte_prob = MonteCarloProblem(prob)
@time sol = solve(monte_prob,dt=0.1,num_monte=10000)

# Spot check the mean and the variance
qs = 0:0.1:1
for i in 2:10
  q = qs[i]
  @test ≈(timestep_mean(sol,i),q,atol=1e-2)
  @test ≈(timestep_meanvar(sol,i)[2],(1-q)*q,atol=1e-2)
end

This test is still stochastic but tends to pass more for 1 and 3 than 2. Its tolerances are set so that way it can go fast enough on Travis and Appveyor, but it should be able to distinguish them more cleanly if you up the number of trajectories.

And of course, you can do this same kind of thing in your own implementation to see. The test could probably be made better if you can get an analytical solution which has a sudden stiffness to induce re-rejections reliably, but at least the Brownian Bridge variation is the most exacting test I have so far.