  The numerical evaluation of statistics plays a crucial role in statistical
physics and its applied fields. It is possible to evaluate the statistics for a
stochastic differential equation with Gaussian white noise via the
corresponding backward Kolmogorov equation. The important notice is that there
is no need to obtain the solution of the backward Kolmogorov equation on the
whole domain; it is enough to evaluate a value of the solution at a certain
point that corresponds to the initial coordinate for the stochastic
differential equation. For this aim, an algorithm based on combinatorics has
recently been developed. In this paper, we discuss a higher-order approximation
of resolvent, and an algorithm based on a second-order approximation is
proposed. The proposed algorithm shows a second-order convergence. Furthermore,
the convergence property of the naive algorithms naturally leads to
extrapolation methods; they work well to calculate a more accurate value with
fewer computational costs. The proposed method is demonstrated with the
Ornstein-Uhlenbeck process and the noisy van der Pol system.
