Convergence of spectra calculated from molecular dynamics trajectories
There are multiple ways in which the VACF, along with other properties derived from molecular dynamics trajectories like the dipole ACF to give IR spectra, can be considered "converged".


The molecular dynamics trajectory itself must be converged by some time-dependent statistical metric, such as measuring fluctuations in the energy or temperature over a number of steps and ensuring the average is below some threshold. This will be related to the ensemble chosen, the thermostat and barostat and associated parameters (if relevant), the integrator used for time evolution, and the size of the simulation, among other things.
Molecular simulations are layers of approximations piled on top of each other. In the case of purely classical dynamics, this usually refers to the force field parameters for the solute and solvent. Force fields are designed for different systems and may or may not be transferable, but the goal is usually to accurately predict dynamics. In the case of the VACF, if the dynamics of a system are described correctly, then the resulting power spectrum should be as well. However, in the case of the dipole ACF, the force field parameters enter directly as the nuclear charges (rather than just indirectly through nuclear velocities), which places stricter requirements on the quality of the force field. One hopes that newer force fields are better, but improvement over time isn't necessarily systematic and there is no quality hierarchy.
As an extension of 2, if the simulation is quantum mechanical in any way (such as in QM/MM, Born-Oppenheimer MD, or Car–Parrinello MD), then there are approximations in how the electronic structure is treated. For spectroscopic properties, it is more clear how to systematically improve results with regard to electronic structure than force fields due to well-established quantum chemical method hierarchies. Benchmarks of this hierarchy have been performed on many properties, including geometric parameters, dipole moments, and vibrational frequencies. In theory, the error in both dynamics and spectra should decrease up going from Hartree-Fock to MP2 to CCSD to CCSD(T), or by using more advanced density functionals (say, BP86 to PBE0 to M06-2X to B2PLYP). In practice, this would require benchmarking to draw any conclusions.
As an extension of 3, there are approximations in QM/MM of how the classically-treated and QM-treated regions interact with each other. Rather than go into more detail, here is a review of different QM/MM methods covering different definitions of the boundary and how one region may interact with the other.


Regarding estimation of the statistical error in resulting spectra, I am only aware of direct comparison (absolute error) between experimental and computed spectra, because it is assumed that the MD trajectory is converged (point 1). Performing multiple independent MD simulations and calculating the VACF from each should give an error measure closer to true statistical error, though I am unsure of how this would be computed.

Here is good almost-review (with good examples) of using Born-Oppenheimer molecular dynamics to calculate power spectra, IR spectra, and Raman spectra. I am not aware of any comparisons between classical MD and BOMD for computing spectra.

On the surface, it may appear as though the introduction of quantum mechanics makes things more difficult; certainly things become more time-consuming, however it is much clearer how results will change by treating electronic structure explicitly and with less severe approximations, rather than with a force field.