Abstract　 　　 

　As a complementary method to N-body simulations, (semi-)analytical models are used to study dark matter subhaloes. Analytical halo merger trees generated using algorithms based on the extended Press-Schechter (EPS) theory provide the backbones of such models. In this work, multiple EPS-based algorithms are tested. The only algorithm in our comparison that yields results indistinguishable with simulations is that by Parkinson et al. (P08). We then develop a new (semi-)analytical model, based on the P08 merger trees and an orbit-averaged subhalo mass-evolution recipe, to study the statistics of dark matter subhaloes. The model accurately reproduces the subhalo mass and velocity functions (SHMFs and SHVFs) from numerical simulations, and has the advantage of being able to probe the halo-to-halo variation of substructure populations at unprecedented precision. We find that the halo-to-halo variation of for Galaxy-sized host haloes, although quite substantial, may not be large enough to salvage the `too big to fail' problem of the Local Group satellites. We also present universal fitting functions for the average SHMF and SHVF that are valid for any host halo mass, at any redshift, and for any LCDM cosmology.