Prosty model liniowy z błędami autokorelowanymi w R [zamknięty]

Jak dopasować model liniowy z błędami autokorelowanymi w R? W stacie użyłbym praispolecenia, ale nie mogę znaleźć odpowiednika R ...

Spójrz na gls(uogólnione najmniejsze kwadraty) z pakietu nlme

Możesz ustawić profil korelacji dla błędów w regresji, np. ARMA, itp:

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

dla błędów ARiMR (1,1).

Oprócz gls()funkcji from nlmemożesz także użyć arima()funkcji w statspakiecie za pomocą MLE. Oto przykład z obiema funkcjami.

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 

Zaletą funkcji arima () jest to, że można dopasować znacznie większą różnorodność procesów błędów ARMA. Jeśli użyjesz funkcji auto.arima () z pakietu prognozy, możesz automatycznie zidentyfikować błąd ARMA:

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

Użyj funkcji gls z pakietu nlme . Oto przykład.

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

Ponieważ model jest montowany z maksymalnym prawdopodobieństwem, musisz podać wartości początkowe. Domyślna wartość początkowa to 0, ale jak zawsze dobrze jest wypróbować kilka wartości, aby zapewnić zbieżność.

Jak zauważył dr G, możesz również użyć innych struktur korelacji, a mianowicie ARMA.

Zauważ, że ogólnie szacunki najmniejszych kwadratów są spójne, jeśli macierz kowariancji błędów regresji nie jest wielokrotnością macierzy tożsamości, więc jeśli dopasujesz model do określonej struktury kowariancji, najpierw musisz sprawdzić, czy jest on odpowiedni.

Możesz użyć przewidywania na wyjściu gls. Widzisz? Prognozuj.gls. Możesz także określić kolejność obserwacji według terminu „forma” w strukturze korelacji. Na przykład:
corr=corAR1(form=~1)wskazuje, że kolejność danych odpowiada kolejności w tabeli. corr=corAR1(form=~Year)wskazuje, że kolejność jest zgodna z współczynnikiem Rok. Wreszcie wartość „0,5” corr=corAR1(0.5,form=~1)?jest na ogół ustawiana na wartość parametru oszacowanego jako reprezentujący strukturę wariancji (phi, w przypadku AR, theta w przypadku MA. .). Opcjonalnie można go skonfigurować i używać do optymalizacji, jak wspomniał Rob Hyndman.