Here is what we/SuanShu suggests:
Determine the lags (p and q) of the ARMA process and fit an ARMA(p, q) model. This is done by the usual ARMA fitting procedure., e.g.,
Select a suitable set of orders (P, Q) for the GARCH process. We can do this by looking at the PACF and ACF of the squared residuals and possibly use Ljung-Box test.
Fit a pure GARCH(P, Q) model to the residuals using conditional MLE.
You can do all steps 1 – 4 in SuanShu by calling the appropriate classes.
See this for more information: