We have developed a skeleton version of a new toolbox for statistical small body dynamics in the Solar system. The propagation parts of the software include perturbations from all major planets, radiation pressure and the Poynting–Robertson effect. Currently, the software is constructed to generate clones of parent bodies taking into account uncertainties in observational parameters and the parent body characteristics. To then sample this distribution in a Monte Carlo fashion. These bodies then release test particles using sublimation models. The parent bodies as well as the particle generation process are described by multivariate probability distributions. In our current usage, the distribution represented orbital elements, critical sublimation radius, density, size and surface activity. The software designed to integrate the released particles over a given time scale and examine close encounters with another body in the solar system. We have examined close encounters with the Earth. We have also created module for calculating orbital similarity functions and to find associations and classifications in data sets. This toolbox is entirely modular enabling the use of every step individually.Validation is performed by simulating known and observed meteor showers, we have simulated the 1933 and 1946 October Draconids as validation, and extended the simulations to the 2011 and 2012 October Draconids. The simulation was performed by ejecting material from comet 21P/Giacobini–Zinner during seven perihelion passages between 1866 and 1972 and propagating the material forward in time. Each perihelion passage was sampled with 50 orbital clones that produced meteoroid streams. In total 850 clones were propagated. The clones were sampled from a multidimensional Gaussian distribution on the orbital elements with width proportional to the given uncertainties. These orbital clones were then sampled from normal distributions on the bulk density, surface activity factor, cometary mass and critical sublimation distance from the Sun, with characteristic values from measurements of 21P/Giacobini–Zinner. Each clone ejected 8,000 particles, each with an individual weight proportional to the mass loss (number of meteoroids) they represented. This generated a total of 6.7 million test particles, out of which 43 thousand entered the Earth's Hill sphere during 1900-2020 and were considered encounters. Using the simulation we produced the unexpected and measured deviation of the meteor mass index from a power low in the 2012 October Draconids, a feature not present in the 2011 October Draconids. We also predict a October Draconids outburst in 2018 with peak on the night between October 8 and 9 that should be larger than the 2011 and 2012 outbursts. Lastly we present some analysis as a proof of concept for the future development of this toolbox.