Multilevel subdivision techniques are presented for the efficient numerical approximation of complicated dynamical behavior. Concretely we develop adaptive methods which allow to extract statistical information on the underlying dynamical system. This is done by an approximation of natural invariant measures as well as (almost) cyclic dynamical components. We discuss issues concerning the implementation (e.g. parallelization strategies) and indicate potential applications of these methods (e.g. to the computation of Lyapunov exponents). The results are illustrated by several numerical examples.