We prove a general theorem on the stochastic convergence of appropriately renormalized models arising from nonlinear stochastic PDEs. The theory of regularity structures gives a fairly automated framework for studying these problems but previous works had to expend significant effort to obtain these stochastic estimates in an ad-hoc manner. In contrast, the main result of this article operates as a black box which automatically produces these estimates for nearly all of the equations that fit within the scope of the theory of regularity structures. Our approach leverages multi-scale analysis strongly reminiscent to that used in constructive field theory, but with several significant twists. These come in particular from the presence of "positive renormalizations" caused by the recentering procedure proper to the theory of regularity structure, from the difference in the action of the group of possible renormalization operations, as well as from the fact that we allow for non-Gaussian driving fields. One rather surprising fact is that although the "canonical lift" is of course typically not continuous on any Hölder-type space containing the noise (which is why renormalization is required in the first place), we show that the "BPHZ lift" where the renormalization constants are computed using the formula given in arXiv:1610.08468, is continuous in law when restricted to a class of stationary random fields with sufficiently many moments.