improved convergence methods are developed for bounded SDEs, it is likely that debiasing will be more competitive with EA. In Giles (2008) it is empir- ically shown in one dimension that using stochastic interpolation techniques with discrete time approximations of SDEs speed up the convergence rates when the SDE is bounded. To our knowledge no theoretical results have followed, and so we did not attempt this here. A number of recent advances enable the use of FKF for a wider class of PDEs than those considered here. In Pollock, Johansen and Roberts (2014) exact algorithms are combined with ϵ -strong simulation in order to estimate boundary crossing times to arbitrary precision. Although this approach is not unbiased, it can be applied to non-polytopal domains. Jenkins (2013) develops EA methodology for reflecting Brownian motion, enabling unbiased estimates of certain Neumann boundary problems. Taylor (2015) develops methodology for exact estimation of local times, extending the potential ap- plication to Neumann problems. Finally, since the original formulation of the Feynman-Kac formulae there have been numerous extensions relating SDEs to other classes of PDEs, for example semilinear PDEs and nonlinear PDEs with forward-backwards SDEs (Pardoux and R˘ a¸scanu, 2014; Pham, 2014), and higher order PDEs with iterated stochastic processes (Thieullen and Vigot, 2015). Simulating from the complementary SDE remains an es- sential step when using these extended formulas, and it is hoped that the work presented here will be built upon in these areas. 27