We study the accuracy of several approximate methods for gravitational dynamics in terms of halo power spectrum multipoles and their estimated covariance matrix. We propagate the differences in covariances into parameter constrains related to growth rate of structure, Alcock-Paczynski distortions and biasing. We consider seven methods in three broad categories: algorithms that solve for halo density evolution deterministically using Lagrangian trajectories (ICE-COLA, Pinocchio and PeakPatch), methods that rely on halo assignment schemes onto dark-matter overdensities calibrated with a target N-body run (Halogen, Patchy) and two standard assumptions about the full density PDF (Gaussian and Lognormal). We benchmark their performance against a set of three hundred N-body simulations, running similar sets of “approximate” simulations with matched initial conditions, for each method. We find that most methods reproduce the monopole to within 5 % 5% , while residuals for the quadrupole are sometimes larger and scale dependent. The variance of the multipoles is typically reproduced within 10 % 10% . Overall, we find that covariances built from ICE-COLA, Pinocchio, PeakPatch, Patchy and the Gaussian approximation yield errors on model parameters within 5 % 5% of those from the N-body based covariance, while for Halogen and lognormal this agreement degrades to ∼ 10 % ∼10% .