GLM with autoregressive term to correct for serial correlation

Question

I have a stationary time series to which I want to fit a linear model with an autoregressive term to correct for serial correlation, i.e. using the formula At = c1*Bt + c2*Ct + ut, where ut = r*ut-1 + et

(ut is an AR(1) term to correct for serial correlation in the error terms)

Does anyone know what to use in R to model this?

Thanks Karl

Answer 1

The GLMMarp package will fit these models. If you just want a linear model with Gaussian errors, you can do it with the arima() function where the covariates are specified via the xreg argument.

Answer 2

There are several ways to do this in R. Here are two examples using the "Seatbelts" time series dataset in the datasets package that comes with R.

The arima() function comes in package:stats that is included with R. The function takes an argument of the form order=c(p, d, q) where you you can specify the order of the auto-regressive, integrated, and the moving average component. In your question, you suggest that you want to create a AR(1) model to correct for first-order autocorrelation in the errors and that's it. We can do that with the following command:

arima(Seatbelts[,"drivers"], order=c(1,0,0),
      xreg=Seatbelts[,c("kms", "PetrolPrice", "law")])

The value for order specifies that we want an AR(1) model. The xreg compontent should be a series of other Xs we want to add as part of a regression. The output looks a little bit like the output of summary.lm() turned on its side.

Another alternative process might be more familiar to the way you've fit regression models is to use gls() in the nlme package. The following code turns the Seatbelt time series object into a dataframe and then extracts and adds a new column (t) that is just a counter in the sorted time series object:

Seatbelts.df <- data.frame(Seatbelts)
Seatbelts.df$t <- 1:(dim(Seatbelts.df)[1])

The two lines above are only getting the data in shape. Since the arima() function is designed for time series, it can read time series objects more easily. To fit the model with nlme you would then run:

library(nlme)
m <- gls(drivers ~ kms + PetrolPrice + law,
         data=Seatbelts.df,
         correlation=corARMA(p=1, q=0, form=~t))
summary(m)

The line that begins with "correlation" is the way you pass in the ARMA correlation structure to GLS. The results won't be exactly the same because arima() uses maximum likelihood to estimate models and gls() uses restricted maximum likelihood by default. If you add method="ML" to the call to gls() you will get identical estimates you got with the ARIMA function above.