I don't know how to select ARMA lag length when doing ARMA-GARCH. Perhaps someone can edit it into this answer.
For the univariate case you want
rugarch package. If you're doing multivariate stuff you want
rmgarch. The reason these are better than other packages is threefold; (i) Support for exogenous variables which I haven't seen in any other package, (ii) support for dynamic conditional correlations, (iii) support for a huge multitude of fGARCH variants.
Let's construct the data to be used as an example. Using $N(0,1)$ will give strange results when you try to use GARCH over it but it's just an example.
data <- rnorm(1000)
We can then compute the ARMA(1,1)-GARCH(1,1) model as an example:
spec <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1),
submodel = NULL, external.regressors = NULL, variance.targeting = FALSE),
mean.model = list(armaOrder = c(1, 1), external.regressors = NULL,
distribution.model = "norm", start.pars = list(), fixed.pars = list()))
garch <- ugarchfit(spec=spec,data=data,solver.control=list(trace=0))
Retrieve ARMA(1,1) and GARCH(1,1) coefficients:
Retrieve time-varying standard deviation:
Retrieve standardized $N(0,1)$ ARMA(1,1) disturbances:
See what else you can pull out of the fit: