Einfaches lineares Modell mit autokorrelierten Fehlern in R [geschlossen]

19

Wie passe ich ein lineares Modell mit autokorrelierten Fehlern in R an? In stata würde ich den praisBefehl verwenden, aber ich kann kein R-Äquivalent finden ...

Zach
quelle

Antworten:

23

Schauen Sie sich gls(generalized least squares) aus dem Paket nlme an

Sie können ein Korrelationsprofil für die Fehler in der Regression festlegen, z. B. ARMA usw.:

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

für ARMA (1,1) -Fehler.

Dr. G
quelle
1
Kann ich die Funktion "Vorhersagen" für einen neuen Datensatz verwenden und dieselbe Fehlerstruktur beibehalten? Woher weiß der Befehl gls, in welcher Reihenfolge sich die Beobachtungen befinden?
Zach
27

Zusätzlich zur gls()Funktion von nlmekönnen Sie die arima()Funktion im statsPaket auch mit MLE verwenden. Hier ist ein Beispiel mit beiden Funktionen.

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 

Der Vorteil der Funktion arima () besteht darin, dass Sie eine viel größere Vielfalt von ARMA-Fehlerprozessen anpassen können. Wenn Sie die Funktion auto.arima () aus dem Prognosepaket verwenden, können Sie den ARMA-Fehler automatisch identifizieren:

require(forecast)    
fit3 <- auto.arima(y, xreg=x)
Rob Hyndman
quelle
1
Was ist das 0.5 für in "corr = corAR1 (0.5, form = ~ 1)"?
Zach
1
Es gibt einen Startwert für die Optimierung. Es macht fast keinen Unterschied, wenn es weggelassen wird.
Rob Hyndman
1
+1 Die arimaOption sieht praisauf den ersten Blick anders aus als die von Stata , ist jedoch flexibler und Sie können sie auch verwenden tsdiag, um ein gutes Bild davon zu erhalten, wie gut Ihre AR (1) -Annahme tatsächlich passt.
Wayne
1
Wie verwende ich arima (), wenn ich nur eine Konstante regressiere?
user43790
7

Verwenden Sie die Funktion gls aus dem Paket nlme . Hier ist das Beispiel.

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

Da das Modell mit maximaler Wahrscheinlichkeit angepasst wird, müssen Sie Startwerte angeben. Der Standard-Startwert ist 0, aber wie immer ist es gut, mehrere Werte auszuprobieren, um die Konvergenz sicherzustellen.

Wie Dr. G ausführte, können Sie auch andere Korrelationsstrukturen verwenden, nämlich ARMA.

Beachten Sie, dass Schätzungen der kleinsten Quadrate im Allgemeinen konsistent sind, wenn die Kovarianzmatrix der Regressionsfehler nicht ein Vielfaches der Identitätsmatrix ist. Wenn Sie also ein Modell mit einer bestimmten Kovarianzstruktur anpassen, müssen Sie zunächst testen, ob dies angemessen ist.

mpiktas
quelle
Was ist das 0.5 für in "corr = corAR1 (0.5, form = ~ 1)"?
Zach
3

Sie können Predict für die gls-Ausgabe verwenden. Siehe "predict.gls". Sie können die Reihenfolge der Beobachtung auch durch den Ausdruck "Form" in der Korrelationsstruktur festlegen. Beispiel:
corr=corAR1(form=~1)Gibt an, dass die Reihenfolge der Daten in der Tabelle angegeben ist. corr=corAR1(form=~Year)gibt an, dass die Reihenfolge der Faktor Jahr ist. Schließlich wird der Wert "0,5" in corr=corAR1(0.5,form=~1)?im Allgemeinen auf den Wert des Parameters gesetzt, der zur Darstellung der Varianzstruktur geschätzt wird (phi im Fall von AR, theta im Fall von MA). .). Es ist optional, es einzurichten und für die Optimierung zu verwenden, wie von Rob Hyndman erwähnt.

Xochitl C.
quelle
Vorhersagen oder angepasste Werte sind sehr unterschiedlich, obwohl sie korrekt sind (gls () und Arima ())? Da gls nur feste Effekte verwendet, bezieht Arima die Arima-Fehler in die angepassten Werte ein.
B_Miner