Sequential Control Variates for Functionals of Markov Processes

Abstract

Using a sequential control variates algorithm, we compute Monte Carlo approximations of solutions of linear partial differential equations connected to linear Markov processes by the Feynman-Kac formula. It includes diffusion processes with or without absorbing/reflecting boundary and jump processes. We prove that the bias and the variance decrease geometrically with the number of steps of our algorithm. Numerical examples show the efficiency of the method on elliptic and parabolic problems. Introduction.We are concerned with the numerical evaluation of E(Ψ(X s : s ≥ t)|X t = x), where (X t ) t is a Markov process (with linear dynamics) and where Ψ belongs to a class of functionals related to Feynman-Kac representations. These issues arise, for example, in physics in the computations of the solution of diffusion equations (see [CDL + 89]), or in finance in the pricing of European options (see [DG95] and the references therein). Monte Carlo methods are usually used to evaluate these expectations for high-dimensional problems or when the functionals are complex. They give a rather poor approximation because of a slow convergence as σ/ √ M , M being the number of simulations and σ 2 the relative variance. A better accuracy can nevertheless be reached by using relevant variance-reduction tools like, for instance, the control variates method or importance sampling [Hal70, New94]. One of the most performing tools is the sequential Monte Carlo approach which consists in using iteratively these variance-reduction ideas [Hal62, Hal70, Boo89]. Using, respectively, importance sampling and control variates, this approach has been recently developed in [BCP00] for Markov chains and in [Mai03] for the numerical integration of multivariate smooth functions. We have introduced in [GM04] a sequential Monte Carlo method to solve the Poisson equation with Dirichlet boundary conditions over square domains. This method was based on Feynman-Kac computations of pointwise solutions combined with a global approximation on Tchebychef polynomials [BM97]. Pointwise solutions were computed using walk-on-spheres (WOS) simulations of stopped Brownian motion, which induces a simulation error due to the absorption layer thickness. We have nevertheless observed a geometric reduction up to a limit of both the simulation error and the variance with the number of steps of the algorithm. The global error was comparable to standard deterministic spectral methods [BM97] while avoiding the resolution of a linear system. Our goal here is twofold:• to extend the scope of the approach to general Markov processes connected to linear elliptic and parabolic Dirichlet problems;

Sequential Monte Carlo methods for diffusion processes

“…Thirdly, it is feasible to extend the ideas here to solve PDEs connected to Markov processes via the Feynman-Kac formula (Gobet & Maire 2005). The SMC ideas, along with different target densities or instrumental distributions, are likely to provide efficient solutions in this context.…”

“…In the context of expectations, the diffusion is Euler discretized and the latter is simulated to yield a Monte Carlo (MC) estimate of the expectation (see Lapeyre et al (2001) for a booklength summary). Various extensions can be found in Gobet & Maire (2005) and Zou & Skeel (2004). When calculating the derivatives, there are two primary methods based on Euler approximation and Mallivian calculus.…”

Some new simulations schemes for the evaluation of Feynman–Kac representations

“…In some situations like spectral methods [4,5] or in the sequential Monte Carlo methods developed in earlier works [13,14,15], the points where the solution is computed are fixed. This means that we can simply replace the simulations by quadrature formulae at some random points of the boundary and of the interior of the domain.…”

“…We describe briefly on the Poisson equation in a domain D an iterative Monte Carlo method introduced in [13,14] which can be used to compute a global approximate solution of many linear elliptic or parabolic partial differential equations. The idea is to use the approximation of the solution at a given step as a control variate for the next step.…”

“…We have introduced sequential Monte Carlo algorithms to compute the solution of the Poisson equation with a great accuracy using either spectral methods [13,14] or domain decomposition methods [15]. In Section 4, we use the formulae of Section 3 to reduce the simulation times and also the number of steps until convergence of these algorithms on two dimensional problems.…”