# 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? – Richard Hardy Feb 24 '17 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 – Darragh Nov 14 '16 at 23:20
• @Darragh, you have to include the future values of time and time^2 via newxreg when forecasting. – Richard Hardy Nov 15 '16 at 6:00
• @Darragh, did it work out? – Richard Hardy Nov 18 '16 at 9:09

Rather than a two stage approach fitting time ,time sq , time cubed there could be level shifts , pulses and untreated ARIMA structure that are giving you "fits" . Following @whuber "Fitting polynomials to data can be a deceptively poor approach: a tiny bit of overfitting can result in models that are grossly bad because higher degree polynomials can (and often do) vary so wildly in between the data values and will be horrible extrapolators". . In building a model that includes memory , level shifts , time trends (as needed .. perhaps a break in trend ! , anomalies , changes in seasonal pulses , parameters and error variance you have all you need under 1 roof , so to speak .