We consider a class of M-estimators of the parameters of a GARCH(p,q) model. These estimators are asymptotically normal, depending on score functions, under milder moment assumptions than the usual quasi maximum likelihood, which makes them more reliable in the presence of heavy tails. We also consider weighted bootstrap approximations of the distributions of these M-estimators and establish their validity. Through extensive simulations, we demonstrate the robustness of these M-estimators under heavy tails and conduct a comparative study of the performance (biases and mean squared errors) of various score functions and the accuracy (confidence interval coverage probabilities) of their bootstrap approximations. In addition to the GARCH(1,1) model, our simulations also involve higher-order models such as GARCH(2,1) and GARCH(1,2) which so far have received relatively little attention in the literature. We also consider the case of order-misspecified models. Finally, we analyze two real financial time series datasets by fitting GARCH(1,1) or GARCH(2,1) models with our M-estimators.