Simple modelo lineal con errores autocorrelated en R

Question

¿Cómo ajustar un modelo lineal con errores autocorrelated en R? En stata utilizaría el prais comando, pero no podemos encontrar un equivalente de R...

Answer 1

Además de la gls() función de nlme, también puede utilizar el arima() función en el stats paquete mediante el MLE. Aquí está un ejemplo con ambas funciones.

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 

La ventaja de la arima() función que pueda adaptarse a una mucho mayor variedad de ARMA de error procesos. Si utilizar el auto.arima() la función de la previsión de paquete, usted puede identificar automáticamente el ARMA de error:

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

Answer 2

Echa un vistazo a gls (mínimos cuadrados generalizados) desde el paquete nlme

Puede establecer un perfil de la correlación de los errores de la regresión, por ejemplo, ARMA, etc.:

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

errores de ARMA(1,1).

Answer 3

Uso de la función gls de paquete nlme. Aquí está el ejemplo.

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

Dado que el modelo está equipado usando máxima verosimilitud usted necesidad de proporcionar a los valores de partida. El valor predeterminado valor inicial es 0, pero como siempre es bueno probar diferentes valores para garantizar la convergencia.

Como el Dr. G señaló también se pueden utilizar otras estructuras de correlación, es decir, ARMA.

Tenga en cuenta que, en general, menos plazas de las estimaciones son consistentes si la matriz de covarianza de los errores de regresión no es múltiplo de matriz de identidad, por lo que si usted ajuste del modelo con los específicos de la estructura de covarianza, en primer lugar usted necesita para comprobar si es adecuada.

Answer 4

Usted puede utilizar predecir en gls de salida. Ver ?predecir.gls. También puede especificar el orden de la observación de la "forma" plazo en la estructura de las correlaciones. Por ejemplo:
corr=corAR1(form=~1) indica que el orden de los datos es el que ellos están en la tabla. corr=corAR1(form=~Year) indica que el orden es el de factor de Año.. Por último, la "0.5" valor en corr=corAR1(0.5,form=~1)? normalmente se configura para el valor del parámetro estimado para representar la variación de la estructura (phi, en el caso de la AR, theta en el caso de MA...). Es opcional para configurar y utilizar para la optimización como Rob Hyndman mencionado.