We consider an unbiased simulation method for multidimensional diffusions based on the parametrix method for solving partial differential equations with Hölder continuous coefficients. This Monte Carlo method which is based on an Euler scheme with random time steps, can be considered as an infinite dimensional extension of the Multilevel Monte Carlo method for solutions of stochastic differential equations with Hölder continuous coefficients. In particular, we study the properties of the variance of the proposed method. In most cases, the method has infinite variance and therefore we propose an importance sampling method to resolve this issue.