The topological heavy-fermion model provides a minimal framework for describing the coexistence of localized moments and itinerant Dirac electrons in magic-angle twisted bilayer graphene (MATBG). Several analytical and numerical methods have been applied to this model; however, whether they provide a realistic description of MATBG remains incompletely understood. In this work, we develop an Hubbard operator approach that incorporates non-local correlations beyond the single-site limit. We benchmark the approximate calculations against numerically exact determinant quantum Monte Carlo simulations of a SU(2)$ symmetric lattice-regularized model. We show that commonly used local approximations, such as Hubbard-I, fail to capture the coupling between localized and itinerant degrees of freedom, leading to incorrect spectral properties in the local-moment regime. In contrast, the Hubbard operator method provides a better description of both correlation functions and spectral features over a regime of parameters, in good agreement with exact numerical methods.