In Monte Carlo estimation of expected values of functionals of solutions of SDEs, we have to generate approximate sample paths, probably by time stepping. Typically, the bias is roughly proportional to a power of the time step, h. The power depends on the time stepper and on the functional. For a functional F(X(T)), so called weak error esimates give first order convergence for the forward Euler and Milstein methods. To bound the bias for more complicated functionals, we study the joint distribution of the path observed at the time step times: (X(h), x(2h), ... X(nh)), where nh=T. We compute the difference between the joint density for the time stepper and the exact joint density. For forward Euler this does not go to zero with h, but for Milstein it goes as sqrt(h). Our analysis leads to Runge Kutta method that is simpler than Milstein but has similar statistical properties. The talk will include background material including the definition of the Milstein method. This is joint work with Peter Glynn and Jose Antonio Perez.