Prediksi 1 langkah lebih maju dengan paket dynlm R

11

Saya sudah cocok dengan model dengan beberapa variabel independen, salah satunya adalah lag dari variabel dependen, menggunakan paket dynlm.

Dengan asumsi saya memiliki ramalan 1 langkah di depan untuk variabel independen saya, bagaimana saya mendapatkan ramalan 1 langkah di depan untuk variabel dependen saya?

Berikut ini sebuah contoh:

library(dynlm)

y<-arima.sim(model=list(ar=c(.9)),n=10) #Create AR(1) dependant variable
A<-rnorm(10) #Create independant variables
B<-rnorm(10)
C<-rnorm(10)
y<-y+.5*A+.2*B-.3*C #Add relationship to independant variables 
data=cbind(y,A,B,C)

#Fit linear model
model<-dynlm(y~A+B+C+L(y,1),data=data)

#Forecast
A<-c(A,rnorm(1)) #Assume we already have 1-step forecasts for A,B,C
B<-c(B,rnorm(1))
C<-c(C,rnorm(1))
y=window(y,end=end(y)+c(1,0),extend=TRUE)
newdata<-cbind(y,A,B,C)
predict(model,newdata)

Dan di sini adalah contoh menggunakan paket dyn, yang berfungsi.

library(dyn)

#Fit linear model
model<-dyn$lm(y~A+B+C+lag(y,-1),data=data)

#Forecast
predict(model,newdata)the dyn packages, which works:
Zach
sumber
Hanya menggunakan dynlmpaket tidak akan memberikan perkiraan untuk variabel dependen Anda. Memberikan perkiraan untuk variabel dependen Anda akan memerlukan model untuk menjelaskannya dan mungkin data tambahan. Saya menyarankan Anda untuk membaca sesuatu tentang regresi multivariat seperti "Analisis Statistik Multivariat Terapan" oleh Johnson dan Wichern. atau kursus perkiraan: duke.edu/~rnau/411home.htm
deps_stats
1
@deps_stats Variabel dependen adalah apa yang ingin saya ramalkan. Saya berasumsi bahwa saya sudah memiliki perkiraan untuk variabel independen saya. Dalam kode contoh saya, y adalah variabel dependen yang saya coba ramalkan, dan A, B, C adalah variabel independen, yang sudah saya perkirakan memiliki ramalan. Jika Anda menjalankan kode contoh yang saya posting, Anda akan memahami sifat masalah saya.
Zach
@Zach: Peringkat Kaggle Bagus! (Saya mengklik tautan di profil mouse-over Anda)
Hugh Perkins

Jawaban:

13

Selamat, Anda telah menemukan bug. Prediksi untuk dynlmdengan data baru rusak jika variabel lagged digunakan. Untuk melihat mengapa melihat output dari

predict(model)
predict(model,newdata=data)

Hasilnya harus sama, tetapi tidak. Tanpa newdataargumen, predictfungsi ini pada dasarnya mengambil modelelemen dari dynlmoutput. Dengan newdataargumen predictmencoba untuk membentuk matriks model baru dari newdata. Karena ini melibatkan parsing formula yang disuplai ke dynlmdan formula memiliki fungsi L, yang hanya didefinisikan fungsi internaly dynlm, matriks model yang salah terbentuk. Jika Anda mencoba men-debug, Anda akan melihat, bahwa variabel dependen tertinggal tidak sedang tertinggal dalam kasus newdataargumen disediakan.

Apa yang dapat Anda lakukan adalah untuk lag variabel dependen dan memasukkannya ke dalam newdata. Berikut adalah kode yang menggambarkan pendekatan ini. Saya menggunakan set.seedsehingga akan mudah direproduksi.

library(dynlm)

set.seed(1)
y<-arima.sim(model=list(ar=c(.9)),n=10) #Create AR(1) dependant variable
A<-rnorm(10) #Create independant variables
B<-rnorm(10)
C<-rnorm(10)
y<-y+.5*A+.2*B-.3*C #Add relationship to independant variables 
data=cbind(y,A,B,C)

#Fit linear model
model<-dynlm(y~A+B+C+L(y,1),data=data)

Inilah perilaku buggy:

> predict(model)
       2        3        4        5        6        7        8        9       10 
3.500667 2.411196 2.627915 2.813815 2.468595 1.733852 2.114553 1.423225 1.470738 
> predict(model,newdata=data)
        1         2         3         4         5         6         7         8         9        10 
2.1628335 3.7063579 2.9781417 2.1374301 3.2582376 1.9534558 1.3670995 2.4547626 0.8448223 1.8762437 

Membentuk newdata

#Forecast fix.
A<-c(A,rnorm(1)) #Assume we already have 1-step forecasts for A,B,C
B<-c(B,rnorm(1))
C<-c(C,rnorm(1))

newdata<-ts(cbind(A,B,C),start=start(y),freq=frequency(y))

newdata<-cbind(lag(y,-1),newdata)
colnames(newdata) <- c("y","A","B","C")

Bandingkan perkiraan dengan model yang cocok:

> predict(model)
       2        3        4        5        6        7        8        9       10 
3.500667 2.411196 2.627915 2.813815 2.468595 1.733852 2.114553 1.423225 1.470738 
> predict(model,newdata=newdata)
       1        2        3        4        5        6        7        8        9       10       11 
      NA 3.500667 2.411196 2.627915 2.813815 2.468595 1.733852 2.114553 1.423225 1.470738 1.102367 

Seperti yang dapat Anda lihat untuk data historis, ramalan tersebut bertepatan dan elemen terakhir berisi ramalan 1 langkah di depan.

mpiktas
sumber
Bagaimana Anda bisa menangani kasus di mana Anda memiliki dua kelambatan dalam rumus yang sama? lag(y,-1)+lag(y,-2)?
Hugh Perkins
1
Nah, maka solusi ini tidak berfungsi. Anda perlu menulis fungsi prediksi Anda sendiri.
mpiktas
Ah, itulah yang saya lakukan sebenarnya:
Hugh Perkins
1
Apakah Anda mempertimbangkan untuk mengirimkannya ke penulis dynlm? Ini adalah situasi yang aneh, bahwa Anda tidak dapat memprediksi menggunakan dynlm.
mpiktas
Hmmm, Anda mengatakan mereka tidak akan secara ajaib memonitor stackoverflow dan memperbaiki bug? Saya kira itu mungkin benar!
Hugh Perkins
2

Mengikuti permintaan @ md-azimul-haque, saya menggali kode sumber saya yang berumur 4 tahun, dan menemukan fungsi bernama tepat berikut ini. Tidak yakin apakah ini yang dicari oleh @ md-azimul-haque?

# pass in training data, test data,
# it will step through one by one
# need to give dependent var name, so that it can make this into a timeseries
predictDyn <- function( model, train, test, dependentvarname ) {
    Ntrain <- nrow(train)
    Ntest <- nrow(test)
    # can't rbind ts's apparently, so convert to numeric first
    train[,dependentvarname] <- as.numeric(train[,dependentvarname])
    test[,dependentvarname] <- NA
    testtraindata <- rbind( train, test )
    testtraindata[,dependentvarname] <- ts( as.numeric( testtraindata[,dependentvarname] ) )
    for( i in 1:Ntest ) {
       cat("predicting i",i,"of",Ntest,"\n")
       result <- predict(model,newdata=testtraindata,subset=1:(Ntrain+i-1))
       testtraindata[Ntrain+i,dependentvarname] <- result[Ntrain + i + 1 - start(result)][1]
    }
    testtraindata <- testtraindata[(Ntrain+1):(Ntrain + Ntest),dependentvarname]
    names(testtraindata) <- 1:Ntest
    return( testtraindata )
}
Hugh Perkins
sumber