Fitting multiple linear regression in R: autocorrelated residuals I'm trying to estimate a multiple linear regression in R with an equation like this:
regr <- lm(rate ~ constant + askings + questions + 0)

askings and questions are quarterly data time-series, constructed with askings <- ts(...).
The problem now is that I got autocorrelated residuals. I know that it is possible to fit the regression using the gls function, but I don't know how to identify the correct AR or ARMA error structure which I have to implement in the gls function. 
I would try to estimate again now with,
gls(rate ~ constant + askings + questions + 0, correlation=corARMA(p=?,q=?))

but I'm unfortunately neither an R expert nor an statistical expert in general to identify p and q.
I would be pleased If someone could give me a useful hint.
Thank you very much in advance!
Jo
 A: Try
library(forecast)
fit <- auto.arima(rate, xreg=cbind(askings,questions))

That will fit the linear model as will as automatically identify an ARMA structure for the errors. It uses MLE rather than GLS, but they are asymptotically equivalent.
A: If prediction is you purpose, you could fit a range of models over parameters:
expand.grid(p = 1:P, q = 1:Q)

where P and Q are the maximal AR(p) and MA(q) terms you wish to include and choose the best fitting model as determined by BIC.
auto.arima() in package forecast will help with this, but it can be coded easily by hand using expand.grid() and loop and the arima() function that comes with R.
The above is fitting on the residuals from a gls() with no correlation structure.
You could also do the whole thing by hand directly with gls() by just fitting lots of models for combinations of p and q and the in built AIC() function.
You could also plot the ACF (acf()) and partial ACF (pacf()) of the residuals from a linear model without correlation structure and use them to suggest the order of model required.
