A Monte Carlo (MC) algorithm is presented for the simulation of the time evolution of aggregation processes featuring multiple components, properties, or conservation laws. Instead of using deterministic differential population balance equations, the MC algorithm utilizes a stochastic approach to aggregation kinetics. As a result, exact simulation of spatially independent aggregation processes is possible without the need for numerical approximations. Furthermore, simulations exactly predict all moments of the size and composition distributions of aggregating particles for both nongelling and gelling kernels and extend these results to the postgelation period. The algorithm is shown to require at most O((Ω 1 Ω 2 …Ω κ ) 1/(κ+1) ) rate-limiting operations per time step for a κ-component aggregation process featuring Ω i monomers of each component i—a substantial performance improvement over the potential of previous methods. Simulation results are presented for bivariate sum, product, and constant kernels, and for the perikinetic (Brownian) kernel.