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 constraints 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 on to dark matter overdensities calibrated with a target N-body run (HALOGEN, PATCHY), and two standard assumptions about the full density probability distribution function (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{{ per cent}}, while residuals for the quadrupole are sometimes larger and scale dependent. The variance of the multipoles is typically reproduced within 10{{ per cent}}. Overall, we find that covariances built from approximate simulations yield errors on model parameters within 10{{ per cent}} of those from the N-body-based covariance.