# How to fit an ARMA process to residuals in R

I have a quadratic regression model with three statistically significant coefficients: intercept, time, and time^2. The residuals of the model show significant autocorrelation. I understand that you can find an ARMA process for the residuals simply by using the ar or arima function and getting the coefficients [i.e. res.ar=ar(resid(fit),method='mle')]... but how do you refit the regression model with the autocorrelated residuals (for purposes in forecasting)? I'm thinking of terms of auto.arima

• I was going through my old answers and noticed this one was not accepted. Do you perhaps need further clarification? Feb 24, 2017 at 14:26

You are correct that it makes sense to estimate the "main part" of the model simultaneously allowing for the error term to have an ARMA structure (as opposed to doing this in two steps, which would be less efficient). Simply add the main regressors (time and time^2) via the argument xreg in the function arima or auto.arima. It will estimate a regression with ARMA errors. You can read more about this and related techniques in an enlightening blog post "The ARIMAX model muddle" by Rob J. Hyndman.
• Ok, so I attempted what you recommended, but I'm still not quite sure if it's correct. I have >adjust=auto.arima(FTEs,xreg=cbind(ftime,I(ftime^2)),max.p=5,max.q=5,trace=T,seasonal=F,allowmean=F) FTEs is the original data. I obtain an ARMA(3,1) process. But I can't seem to make any forecasts with the model. Instead, I get an error: Error in predict.Arima(adjust) : 'xreg' and 'newxreg' have different numbers of columns Nov 14, 2016 at 23:20
• @Darragh, you have to include the future values of time and time^2 via newxreg when forecasting. Nov 15, 2016 at 6:00