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...
 A: 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.
A: 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)

A: You can use predict on gls output. See ?predict.gls. Also you can specify the order of the observation by the "form" term in the correlation structure. 
For example:
corr=corAR1(form=~1) indicates that order of the data is the one they are in the table.
corr=corAR1(form=~Year) indicates that the order is the one of factor Year..
Finally the "0.5" value in corr=corAR1(0.5,form=~1)? is generally set up to the value of the parameter estimated to represent the variance structure (phi, in case of AR, theta in case of MA...). It is optional to set it up and use for optimization as Rob Hyndman mentioned.
A: 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.
