2009-12-14 7 views
7

Ho una serie temporale stazionaria a cui voglio adattare un modello lineare con un termine autoregressivo per correggere la correlazione seriale, cioè utilizzando la formula At = c1 * Bt + c2 * Ct + ut, dove ut = r * ut-1 + etGLM con termine autoregressivo per correggere la correlazione seriale

(ut è un AR (1) termine per correggere autocorrelazione nei termini di errore)

qualcuno sa cosa usare in R per modellare Questo?

Grazie Karl

risposta

8

Il pacchetto GLMMarp si adatta a questi modelli. Se si desidera un modello lineare con errori gaussiani, è possibile farlo con la funzione arima() in cui le covariate vengono specificate tramite l'argomento xreg.

+0

Il pacchetto GLMMarp è stato rimosso dal repository CRAN. Conoscete qualche altro pacchetto che porterebbe a termine questo? – TMS

+0

Tutti i pacchetti che conosco sono elencati su http://cran.r-project.org/web/views/TimeSeries.html –

2

Qual è la sua funzione di collegamento?

Il modo in cui lo descrivi sembra una regressione lineare di base con errori autocorrelati. In tal caso, un'opzione è utilizzare lm per ottenere una stima coerente dei coefficienti e utilizzare Newey-West HAC standard errors.

Non sono certo la risposta migliore per GLM in generale.

3

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.