In this work, we analyze the Born, Bogoliubov, Green, Kirkwood, and Yvon (BBGKY) hierarchy of equations for describing the full time evolution of a many-body fermionic system in terms of its reduced density matrices (at all orders). We provide an exhaustive study of the challenges and open problems linked to the truncation of such a hierarchy of equations to make them practically applicable. We restrict our analysis to the coupled evolution of the one- and two-body reduced density matrices, where higher-order correlation effects are embodied into the approximation used to close the equations. We prove that within this approach, the number of electrons and total energy are conserved, regardless of the employed approximation. Further, we demonstrate that although most of the truncation schemes available in the literature give acceptable ground-state energy, when applied to describe driven electron dynamics, they exhibit undesirable and unphysical behavior, e.g., violation and even divergence in local electronic density, both in weakly and strongly correlated regimes. We illustrate and analyze these problems within the few-site Hubbard model. The model can be solved exactly and provides a unique reference for our detailed study of electron dynamics for different values of interaction, different initial conditions, and the large set of approximations considered here. Moreover, we study the role of compatibility between two hierarchical equations and positive semidefiniteness of reduced density matrices in the instability of electron dynamics. We show that even if the used approximation holds the compatibility, electron dynamics can still diverge when positive definitiveness is violated. We propose some partial solutions of such problems and point out the main paths for future work in order to make this approach applicable for the description of the correlated electron dynamics in complex systems.