# Simple linear model with autocorrelated errors in R

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...

-

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.

-
Can I use the "predict" function on a new dataset, and retain the same error structure? How does the gls command know what order the observations are in? –  Zach Jan 24 '11 at 0:23

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)

-
What is the 0.5 for in "corr=corAR1(0.5,form=~1)?" –  Zach Jan 24 '11 at 0:30
It gives a starting value for the optimization. It makes almost no difference if it is omitted. –  Rob Hyndman Jan 24 '11 at 2:11
+1 The arima option looks more different from Stata's prais at first glance, but it's more flexible and you can also use tsdiag to get a nice visual of how well your AR(1) assumption actually fits. –  Wayne Aug 24 at 14:37

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.

-
What is the 0.5 for in "corr=corAR1(0.5,form=~1)?" –  Zach Jan 24 '11 at 0:29