Currently you have a two-step estimation procedure: first fit the linear model, then fit a model for the residuals. A more efficient approach would be to estimate the model jointly. Then you would be able to forecast easier and hypothesis testing would also be more straightforward because appropriate standard errors of the linear model would be readily available (unlike in the two-step approach).
Joint estimation can be done via the
forecast::auto.arima functions in R. The former would require you to specify the error structure manually, while the latter would pick the error structure according to some information criterion (AICc by default; alternatively AIC or BIC). You would specify the linear model by including the regressors via the optional argument
y is the dependent variable and
X is your design matrix (the matrix of regressors).
Forecasting from jointly estimated model
Once you have estimated your model, you could apply the
predict method onto it to get forecasts. You would have to supply new values of regressors via argument
newxreg, while the errors would be forecasted automatically given the estimated model object. E.g.
Xnew are the new values for the regressors.
Forecasting from model estimated in two steps
If you are to forecast from the separately estimated linear model and the error equation, here is how it goes:
- Forecast for the linear model:
Obtain new values of your independent variables and multiply them to the estimated coefficient from the linear model (
[-1] deal with the intercept if it is not included within
- Forecast for the error equation:
predict function for the estimated MA (or ARMA) model (
- Add the two forecasts together: