The two-loop anomalous dimension of the chiral matter superfields is calculated for a general N = 1 supersymmetric gauge theory regularized by higher covariant derivatives. We obtain both the anomalous dimension defined in terms of the bare couplings, and the one defined in terms of the renormalized couplings for an arbitrary renormalization prescription. For the one-loop finite theories we find a simple relation between the higher derivative regulators under which the anomalous dimension defined in terms of the bare couplings vanishes in the considered approximation. In this case the one-loop finite theory is also two-loop finite in the HD+MSL scheme. Using the assumption that with the higher covariant derivative regularization the NSVZ equation is satisfied for RGFs defined in terms of the bare couplings, we construct the expression for the three-loop β-function. Again, the result is written both for the β-function defined in terms of the bare couplings and for the one defined in terms of the renormalized couplings for an arbitrary renormalization prescription.