Prédictions à 1 pas avec le package dynlm R

11

J'ai ajusté un modèle avec plusieurs variables indépendantes, dont l'une est le décalage de la variable dépendante, en utilisant le package dynlm.

En supposant que j'ai des prévisions à un pas pour mes variables indépendantes, comment puis-je obtenir des prévisions à un pas pour mes variables dépendantes?

Voici un exemple:

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)

Et voici un exemple d'utilisation du package dyn, qui fonctionne.

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
la source
Utiliser uniquement le dynlmpackage ne fournira pas de prévisions pour vos variables dépendantes. Fournir des prévisions pour vos variables dépendantes nécessitera un modèle pour les expliquer et probablement des données supplémentaires. Je vous suggère de lire quelque chose sur la régression multivariée comme "Applied Multivariate Statistical Analysis" de Johnson et Wichern. ou un cours sur les prévisions: duke.edu/~rnau/411home.htm
deps_stats
1
@deps_stats La variable dépendante est ce que je veux prévoir. Je suppose que j'ai déjà des prévisions pour mes variables indépendantes. Dans mon exemple de code, y est la variable dépendante que j'essaie de prévoir, et A, B, C sont les variables indépendantes pour lesquelles j'ai déjà des prévisions. Si vous exécutez l'exemple de code que j'ai publié, vous comprendrez la nature de mon problème.
Zach
@Zach: Belle cote Kaggle! (J'ai cliqué sur le lien dans votre profil de survol)
Hugh Perkins

Réponses:

13

Félicitations, vous avez trouvé un bug. La prédiction pour les dynlmnouvelles données est rompue si des variables retardées sont utilisées. Pour voir pourquoi regarder la sortie de

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

Les résultats devraient être les mêmes, mais ils ne le sont pas. Sans newdataargument, la predictfonction récupère essentiellement l' modelélément de la dynlmsortie. Avec newdataargument predictessaie de former une nouvelle matrice de modèle à partir de newdata. Étant donné que cela implique l'analyse de la formule fournie à dynlmet que la formule a une fonction L, qui n'est définie qu'en interne dynlm, la matrice de modèle incorrecte est formée. Si vous essayez de déboguer, vous verrez que la variable dépendante retardée n'est pas retardée dans le cas où l' newdataargument est fourni.

Ce que vous pouvez faire est de retarder la variable dépendante et de l'inclure dans le newdata. Voici le code illustrant cette approche. J'utilise set.seeddonc ce serait facilement reproductible.

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)

Voici le comportement du 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 

Former le 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")

Comparer les prévisions avec l'ajustement du modèle:

> 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 

Comme vous pouvez le voir pour les données historiques, les prévisions coïncident et le dernier élément contient les prévisions à 1 étape.

mpiktas
la source
Comment pouvez-vous gérer le cas où vous avez deux retards dans la même formule? lag(y,-1)+lag(y,-2)?
Hugh Perkins
1
Eh bien, cette solution ne fonctionne pas. Vous devez écrire votre propre fonction de prédiction.
mpiktas
Ah, c'est ce que j'ai fait en fait :-P
Hugh Perkins
1
Avez-vous pensé à le soumettre aux auteurs de dynlm? C'est une situation bizarre, que vous ne pouvez pas prévoir en utilisant dynlm.
mpiktas
Hmmm, vous dites qu'ils ne vont pas surveiller par magie le débordement de pile et corriger les bugs? Je suppose que c'est probablement vrai!
Hugh Perkins
2

Suite à la demande de @ md-azimul-haque, j'ai fouillé dans mon code source de 4 ans et j'ai trouvé la fonction suivante nommée de manière appropriée. Vous ne savez pas si c'est ce que recherche @ 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
la source