ARMA/GARCH estimation in sequence I have a time series that shows a nonstationary seasonal autoregressive component as well as known heteroshedasticity. In order to model the series, I have fit a seasonal ARIMA model for the mean with the auto.arima model in the forecast package in R and a GARCH model on the residuals of the ARIMA model. 
Is the procedure of sequentially estimating ARIMA and GARCH model correct or would it have been better to jointly model the mean and the variance of the series? If this were correct, is there a (possibly R) function to do it?
 A: Doing joint estimation is the preferred way. If you do estimation in two stages, a logical inconsistency arises. In the first stage of seasonal ARIMA estimation there is an effective assumption of conditional homoskedasticity. It is contradicted in the second stage when you explicitly model conditional heteroskedasticity using a GARCH model.
If you have no MA terms in the ARIMA model, you will still get consistent parameter estimates even when neglecting GARCH errors, but these estimates will be inefficient. If you do have MA terms in the ARIMA model, the estimates of the parameters in the ARIMA model will not even be consistent.
Functions ugarchspec and ugarchfit in package rugarch (see here for a vignette)  allow specifying and estimating ARMA+GARCH models simultaneously for a variety of GARCH model classes. Unfortunately, seasonal ARMA models do not seem to be implemented there. Perhaps you could try seasonally adjusting your series before fitting an ARMA+GARCH model (although this would be suboptimal if the "true" model is seasonal ARIMA with conditionally heteroskedastic errors).
