We use stellar dynamics as a testbed for statistical closure theory. We focus on the process of « Vector Resonant Relaxation, » a long-range, non-linear, and correlated relaxation mechanism that drives the reorientation of stellar orbital planes around a supermassive black hole. This process provides a natural setting to evaluate the predictive power of generic statistical closure schemes for dynamical correlation functions, in the fully non-linear and non-perturbative regime. We develop a numerical scheme that explicitly implements the seminal « Martin-Siggia-Rose » formalism at one-loop order via an iterative fixed-point approach, thereby improving upon the bare order from the « Direct Interaction Approximation. » Using this framework, we quantitatively validate the ability of the formalism to predict (i) the two-point two-time correlation function; (ii) the renormalised three-point interaction vertex; (iii) the three-point three-time correlation function. These predictions are compared to direct measurements from numerical simulations. We conclude by discussing the limitations of this approach and presenting possible future venues.