Modelo linear simples com erros autocorrelacionados em R [fechado]

19

Como encaixo um modelo linear com erros autocorrelacionados em R? Em stata, eu usaria o praiscomando, mas não consigo encontrar um equivalente em R ...

Zach
fonte

Respostas:

23

Veja gls(mínimos quadrados generalizados) do pacote nlme

Você pode definir um perfil de correlação para os erros na regressão, por exemplo, ARMA, etc:

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

para erros ARMA (1,1).

Dr G
fonte
1
Posso usar a função "prever" em um novo conjunto de dados e manter a mesma estrutura de erro? Como o comando gls sabe em que ordem as observações estão?
Zach
27

Além da gls()função de nlme, você também pode usar a arima()função no statspacote usando o MLE. Aqui está um exemplo com as duas funções.

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 

A vantagem da função arima () é que você pode ajustar uma variedade muito maior de processos de erro do ARMA. Se você usar a função auto.arima () do pacote de previsão, poderá identificar automaticamente o erro ARMA:

require(forecast)    
fit3 <- auto.arima(y, xreg=x)
Rob Hyndman
fonte
1
Qual é o valor 0,5 em "corr = corAR1 (0,5, formulário = ~ 1)?"
Zach
1
Ele fornece um valor inicial para a otimização. Não faz quase nenhuma diferença se for omitido.
Rob Hyndman
1
+1 A arimaopção parece mais diferente da Stata praisà primeira vista, mas é mais flexível e você também pode usar tsdiagpara obter um bom visual de quão bem sua suposição de AR (1) realmente se encaixa.
Wayne
1
Como uso arima () se regredir y apenas em uma constante?
precisa saber é o seguinte
7

Use a função gls do pacote nlme . Aqui está o exemplo.

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

Como o modelo é ajustado usando a probabilidade máxima, você precisa fornecer valores iniciais. O valor inicial padrão é 0, mas como sempre, é bom tentar vários valores para garantir a convergência.

Como o Dr. G apontou, você também pode usar outras estruturas de correlação, a saber, o ARMA.

Observe que, em geral, as estimativas de mínimos quadrados são consistentes se a matriz de covariância dos erros de regressão não for múltipla da matriz de identidade; portanto, se você ajustar o modelo a uma estrutura de covariância específica, primeiro precisará testar se é apropriado.

mpiktas
fonte
Qual é o valor 0,5 em "corr = corAR1 (0,5, formulário = ~ 1)?"
Zach
3

Você pode usar prever na saída GLS. Veja? Predict.gls. Também é possível especificar a ordem da observação pelo termo "formulário" na estrutura de correlação. Por exemplo:
corr=corAR1(form=~1)indica que a ordem dos dados é a que eles estão na tabela. corr=corAR1(form=~Year)indica que a ordem é a do fator Ano. Finalmente, o valor "0,5" corr=corAR1(0.5,form=~1)?é geralmente definido como o valor do parâmetro estimado para representar a estrutura de variação (phi, no caso de AR, teta, no caso de MA .. .). É opcional configurá-lo e usá-lo para otimização, como Rob Hyndman mencionou.

Xochitl C.
fonte
Os valores preditivos ou ajustados serão muito diferentes, apesar de corretos (gls () versus Arima ())? Como o gls usará apenas efeitos fixos, enquanto o Arima incluirá os erros de arima nos valores ajustados.
B_Miner 4/15