We present a self-consistent multilevel simulation method to estimate equilibrium properties of multiscale physical systems. Inspired by multilevel Markov chain Monte Carlo methods, we use an efficient but inexact coarse-grained approximation as the starting point of a hierarchical method to simulate the exact system of interest. We apply this method to highly size-asymmetric binary mixtures of hard spheres, where the scale separation between big and small particles poses a substantial challenge to standard Monte Carlo simulation methods. The big particles of this system display interesting collective behaviour, including a de-mixing transition at large size-ratios and high small-particle densities. To investigate this system, we first develop a two-level method that enables us to simulate the system up to this transition, providing the first computational evidence of its existence and locating the associated critical point. Subsequently, we discuss a generalisation of the two-level method that makes use of multiple levels of physical coarse-graining, and apply a three-level version to the binary mixture. For this example, we compare the numerical and asymptotic performance of the two- and three-level method. We show that taking an intermediate level into account can lower the variance of the method at fixed computational cost.