Applying GARCH Model to Forecasts I'm reading up on GARCH models in Springer Introductory Time Series, and had a question on how we actually apply the model to forecasts.
The data can be grabbed like so:
stemp<- scan("http://www.maths.adelaide.edu.au/andrew.metcalfe/Data/stemp.dat")
stemp.ts<- ts(stemp, start = 1850, frequency = 12)

I originally fit a seasonal ARIMA to the data using auto.arima. The squared residuals appear to be correlated so I thought I'd fit a GARCH(1, 1) model to the residuals.
library(forecast)
stemp.arima<- auto.arima(stemp.ts)
acf(stemp.arima$residuals^2)

library(tseries)
stemp.garch<- garch(stemp.arima$residuals, trace = F)

The correlogram of the squared stemp.garch residuals appear to be white noise. Now that I have that knowledge, though, I don't understand how it would impact my predictions. I was thinking I could simulate a GARCH process using the coefficients from the GARCH output and the distribution of the seasonal ARIMA residuals. I could then take that and add it to my seasonal ARIMA predictions. Not sure if that's a good method, though.
 A: To gain efficiency, it would make sense to fit a SARIMA-GARCH model in one step rather than starting with SARIMA and following it up with GARCH. Unfortunately, I am not aware of any R package that would have the functionality for one-step fitting of SARIMA-GARCH models. Restricted ARIMA-GARCH could be used as a substitute, though.
In the case of a two-step solution, having a volatility forecast from GARCH allows you to replace the constant $h$-step-ahead forecast intervals with time-varying ones. You would derive them from the predicted conditional variance from the GARCH model. 


*

*For example, the one-step-ahead 80% forecast interval would be the point forecast from SARIMA plus the interval between 10% and 90% quantiles of the forecasted error distribution from GARCH. 

*For $h$-step-ahead forecast intervals when $h>1$, you would estimate the variance as the sum of conditional variances over the $h$ periods and use that to determine the expected forecast error distribution from GARCH. Then take the point forecast and add quantiles as in the case of $h=1$.

