Ci sono diversi modi per fare questo in R. Ecco due esempi che utilizzano il "Seatbelts" time series dataset in the datasets package che viene fornito con R.
La funzione arima()
è disponibile in package: le statistiche che è incluso con R. La funzione prende un argomento di il modulo order=c(p, d, q)
in cui è possibile specificare l'ordine del componente auto-regressivo, integrato e della media mobile. Nella tua domanda, suggerisci di voler creare un modello AR (1) per correggere l'autocorrelazione del primo ordine negli errori e il gioco è fatto. Possiamo farlo con il seguente comando:
arima(Seatbelts[,"drivers"], order=c(1,0,0),
xreg=Seatbelts[,c("kms", "PetrolPrice", "law")])
Il valore per ordine specifica che vogliamo un modello AR (1). Il componente xreg dovrebbe essere una serie di altre X che vogliamo aggiungere come parte di una regressione. L'output assomiglia un po 'all'uscita di summary.lm()
ruotata su un lato.
Un altro processo alternativo potrebbe essere più familiare al modo in cui i modelli di regressione si sono adattati utilizzando gls()
nello nlme package. Il codice seguente trasforma l'oggetto serie temporale delle cinture in dataframe e poi estratti e aggiunge una nuova colonna (t) che è solo un contatore nell'oggetto serie temporale ordinato:
Seatbelts.df <- data.frame(Seatbelts)
Seatbelts.df$t <- 1:(dim(Seatbelts.df)[1])
Le due righe sopra sono solo ottenere i dati in forma. Poiché la funzione arima()
è progettata per le serie temporali, può leggere più facilmente gli oggetti delle serie temporali. Per adattare il modello con NLME si sarebbe poi eseguire:
library(nlme)
m <- gls(drivers ~ kms + PetrolPrice + law,
data=Seatbelts.df,
correlation=corARMA(p=1, q=0, form=~t))
summary(m)
La riga che inizia con "correlazione" è il modo in cui si passa nella struttura di correlazione ARMA a GLS. I risultati non saranno esattamente gli stessi perché arima()
utilizza la massima verosimiglianza per stimare i modelli e gls()
utilizza la massima verosimiglianza per impostazione predefinita. Se aggiungi method="ML"
alla chiamata a gls()
, otterrai le stesse stime ottenute con la funzione ARIMA sopra.
Il pacchetto GLMMarp è stato rimosso dal repository CRAN. Conoscete qualche altro pacchetto che porterebbe a termine questo? – TMS
Tutti i pacchetti che conosco sono elencati su http://cran.r-project.org/web/views/TimeSeries.html –