The accurate simulation of time-dependent many-body systems is among the most challenging topics in modern theoretical physics. For many-body systems the number of degrees of freedom is large and conservation laws are rare making the numerical solution very demanding. Quantum mechanics adds further complications. It requires the collective description of the system in terms of the many-body wavefunction whose time evolution is governed by the Schrödinger equation. The collective nature of the description forbids to decompose the system into individual entities and makes the calculation of a quantum many-body system particularly complicated. In fact, the direct solution of the Schrödinger equation is not a viable option for atomic and molecular systems consisting of more than two electrons. To overcome this limitation a variety of approaches have been developed for the calculation of ground state properties of multi-electron systems relevant for the kinetics of chemical reactions. These methods of quantum chemistry can be roughly divided into two groups: methods that are based on the wavefunction as the fundamental object, such as multiconfigurational Hartree-Fock, and methods based on reduced quantities, such as density functional theory. While wavefunction-based methods are very accurate their applicability is limited to small systems due to the exponential scaling with particle number. Density functional theory, on the other hand, can treat large and extended systems, however, at the price of introducing the exchange-correlation functional whose exact form is unknown and whose approximations are hard to improve systematically. Time-dependent formulations have been achieved for both of the aforementioned approaches. The multiconfigurational time-dependent Hartree-Fock method is among the most accurate approaches to simulate dynamical many-body systems. Time-dependent density functional theory is regularly employed to simulate time-dependent large-scale systems. In this thesis, we aim for a method that combines the best of both worlds, the accuracy of wavefunctionbased approaches with the efficiency of polynomial scaling as in density functional theory. To this end we propagate the time-dependent two-particle reduced density matrix. As a hybrid between the electron density and the many-body wavefunction the two-particle reduced density matrix fully includes two-particle correlations which is a prerequisite to accurately capture effects that arise from electron-electron interactions. Further observables such as kinetic energy spectra, ionization probabilities or the total energy can be expressed directly without invoking approximate read-out functionals as required in density functional theory. We develop a closed equation of motion for the two-particle reduced density matrix by constructing a novel reconstruction functional for the three-particle reduced density matrix that preserves norm, energy, and spin symmetries during time propagation. Further, we show how to avoid instabilities associated with the violation of N-representability that have been a considerable limitation in previous approaches. We have implemented the time-dependent two-particle reduced density matrix method to describe high-harmonic generation from fully three-dimensional multi-electron atoms as well as from one-dimensional molecules. We benchmark the performance of the time-dependent two-particle reduced density matrix method by comparing it to a state of the art multiconfigurational time-dependent Hartree-Fock calculation as well as to time-dependent density functional theory. We find very good agreement between the time-dependent two-particle reduced density matrix method and the multiconfigurational time-dependent Hartree-Fock method while time-dependent density functional theory within the local density approximation shows clear deviations indicating that the correct treatment of two-particle correlations is essential to obtain accurate high-harmonic spectra.