<p>This study presents the stochastic multibody dynamics formulation to investigate the long time evaluation of dynamic systems under uncertainty. The generalized polynomial chaos method has been developed to investigate uncertainty in the fields of structural dynamics and fluid dynamics. However, in multibody dynamics fields, uncertainty analysis is challenging due to the high degree of nonlinearity of multibody systems and the typical problem of polynomial chaos, which degrades accuracy in the time domain. To address these challenges, a sequential orthogonalization technique with whitening transformation is employed to resolve the ill-condition issue of the Gram-Schmidt process. A correlation based selection of model solutions is also considered to improve computational efficiency at each updating step as well as numerical stability. These methodologies improve polynomial basis reconstruction and manage statistical dependencies, thereby mitigating the rapid accuracy decay over time and ensuring high precision for extended duration. Numerical experiments demonstrate that the proposed method significantly out-performs conventional polynomial chaos methods.</p>