Simple linear model with autocorrelated errors in R

Question

How do I fit a linear model with autocorrelated errors in R? In stata I would use the prais command, but I can't find an R equivalent...

Answer 1

Have a look at gls (generalized least squares) from the package nlme

You can set a correlation profile for the errors in the regression, e.g. ARMA, etc:

 gls(Y ~ X, correlation=corARMA(p=1,q=1))

for ARMA(1,1) errors.

Answer 2

In addition to the gls() function from nlme, you can also use the arima() function in the stats package using MLE. Here is an example with both functions.

x <- 1:100
e <- 25*arima.sim(model=list(ar=0.3),n=100)
y <- 1 + 2*x + e

###Fit the model using gls()
require(nlme)
(fit1 <- gls(y~x, corr=corAR1(0.5,form=~1)))
Generalized least squares fit by REML
  Model: y ~ x 
  Data: NULL 
  Log-restricted-likelihood: -443.6371

Coefficients:
(Intercept)           x 
   4.379304    1.957357 

Correlation Structure: AR(1)
 Formula: ~1 
 Parameter estimate(s):
      Phi 
0.3637263 
Degrees of freedom: 100 total; 98 residual
Residual standard error: 22.32908 

###Fit the model using arima()
(fit2 <- arima(y, xreg=x, order=c(1,0,0)))

Call:
arima(x = y, order = c(1, 0, 0), xreg = x)

Coefficients:
         ar1  intercept       x
      0.3352     4.5052  1.9548
s.e.  0.0960     6.1743  0.1060

sigma^2 estimated as 423.7:  log likelihood = -444.4,  aic = 896.81 

The advantage of the arima() function is that you can fit a much larger variety of ARMA error processes. If you use the auto.arima() function from the forecast package, you can automatically identify the ARMA error:

require(forecast)    
fit3 <- auto.arima(y, xreg=x)

Answer 3

Use function gls from package nlme. Here is the example.

##Generate data frame with regressor and AR(1) error. The error term is 
## \eps_t=0.3*\eps_{t-1}+v_t
df <- data.frame(x1=rnorm(100), err=filter(rnorm(100)/5,filter=0.3,method="recursive"))

##Create ther response
df$y <- 1 + 2*df$x + df$err

###Fit the model
gls(y~x, data=df, corr=corAR1(0.5,form=~1))

Generalized least squares fit by REML
  Model: y ~ x 
  Data: df 
  Log-restricted-likelihood: 9.986475

 Coefficients:
 (Intercept)           x 
   1.040129    2.001884 

 Correlation Structure: AR(1)
 Formula: ~1 
 Parameter estimate(s):
     Phi 
 0.2686271 
Degrees of freedom: 100 total; 98 residual
Residual standard error: 0.2172698 

Since model is fitted using maximum likelihood you need to supply starting values. The default starting value is 0, but as always it is good to try several values to ensure the convergence.

As Dr. G pointed out you can also use other correlation structures, namely ARMA.

Note that in general least squares estimates are consistent if covariance matrix of regression errors is not multiple of identity matrix, so if you fit model with specific covariance structure, first you need to test whether it is appropriate.