7
votes
$\begingroup$

I have a dataset that I want to fit a simple linear model to, but I want to include the lag of the dependent variable as one of the regressors. Then I want to predict future values of this time series using forecasts I already have for the independent variables. The catch is: how do I incorporate the lag into my forecast?

Here's an example:

#A function to calculate lags
lagmatrix <- function(x,max.lag){embed(c(rep(NA,max.lag),x),max.lag)}
lag <- function(x,lag) {
 out<-lagmatrix(x,lag+1)[,lag]
 return(out[1:length(out)-1])
}

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

#Fit linear model
lag1<-lag(y,1)
model<-lm(y~A+B+C+lag1)
summary(model)

#Forecast linear model
A<-rnorm(50) #Assume we know 50 future values of A, B, C
B<-rnorm(50)
C<-rnorm(50)
lag1<-  #################This is where I'm stuck##################

newdata<-as.data.frame(cbind(A,B,C,lag1))
predict.lm(model,newdata=newdata)
$\endgroup$
2
votes
$\begingroup$

It's clear the solution I posted previously is inadequate and inelegant. Here is my second attempt, which 100% solves my problem. Please let me know if you spot any bugs! I will cross post to stack overflow, if you all think that would be a better place to get comments on my code.

#A function to iteratively predict a time series
ipredict <-function(model, newdata, interval = "none",
        level = 0.95, na.action = na.pass, weights = 1) {
    P<-predict(model,newdata=newdata,interval=interval,
        level=level,na.action=na.action,weights=weights)
    for (i in seq(1,dim(newdata)[1])) {
        if (is.na(newdata[i])) {
            if (interval=="none") {
                P[i]<-predict(model,newdata=newdata,interval=interval,
                    level=level,na.action=na.action,weights=weights)[i]
                newdata[i]<-P[i]
            }
            else{
                P[i,]<-predict(model,newdata=newdata,interval=interval,
                    level=level,na.action=na.action,weights=weights)[i,]
                newdata[i]<-P[i,1]
            }
        }
    }
    P_end<-end(P)[1]*frequency(P)+(end(P)[2]-1) #Convert (time,period) to decimal time
    P<-window(P,end=P_end-1*frequency(P)) #Drop last observation, which is NA
    return(P)
}


#Example usage:
library(dyn)
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)
Error<-rnorm(10)
y<-y+.5*A+.2*B-.3*C+.1*Error #Add relationship to independant variables 
data=cbind(y,A,B,C)

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

#Forecast linear model
A<-c(A,rnorm(5))
B<-c(B,rnorm(5))
C<-c(C,rnorm(5))
y=window(y,end=end(y)+c(5,0),extend=TRUE)
newdata<-cbind(y,A,B,C)
P1<-ipredict(model.dyn,newdata)
P2<-ipredict(model.dyn,newdata,interval="prediction")

#Plot
plot(y)
lines(P1,col=2)
$\endgroup$
  • $\begingroup$ @Zach, did you look at the results? Your function still does not predict the future values. In the loop you are basicaly repeating the same thing you done before the loop. Also note, that your function will fail if the dependent variable is not in the first column of newdata. $\endgroup$ – mpiktas Feb 1 '11 at 8:22
  • $\begingroup$ @mpiktas: what I am trying to do is this: 1. Starting with the first NA value, use all previous data to construct a 1-step forecast using dyn. Replace the NA value with this forecast, and repeat. I thought this would iteratively step through all future NA values... what am I missing? $\endgroup$ – Zach Feb 1 '11 at 14:21
  • $\begingroup$ @Zach, the problem is with one-step forecast. Apparently it does not work. $\endgroup$ – mpiktas Feb 1 '11 at 14:30
  • $\begingroup$ @It seems to work with the dyn function, rather than the dynlm function. $\endgroup$ – Zach Feb 1 '11 at 16:17
  • $\begingroup$ @Zach, I thought dyn$lm was a typo of dynlm. I did not know about package dyn. With it everything works. For the code, my other comment about first column of newdata still holds. Also you might want to look into using call objects. They make the passing of default arguments much more easier. You can ask about it in stackoverflow. $\endgroup$ – mpiktas Feb 2 '11 at 7:52
2
votes
$\begingroup$

One more method, which has been suggested in other topics, is to just use the arima function with xregs. Arima seems to be able to make forecasts from a new set of xregs just fine.

$\endgroup$
0
votes
$\begingroup$

Ok, I answered my own problem, but my solution could use more testing and probably isn't perfect. Suggestions would be appreciated!

First of all I used a modified version of the parseCall function available here:

parseCall <- function(obj) {
    if (class(obj) != 'call') {
        stop("Must supply a 'call' object")
    }

    srep <- deparse(obj)
    if (length(srep) >1) srep <- paste(srep,sep='',collapse='')

    fname <- unlist(strsplit(srep,"\\("))[1]

    func <- unlist(strsplit(srep, paste(fname,"\\(",sep='')))[2]
    func <- unlist(strsplit(func,""))
    func <- paste(func[-length(func)],sep='',collapse="")

    func <- unlist(strsplit(func,","))

    vals <- list()
    nms <- c()
    cnt <- 1
    for (args in func) {
        arg <- unlist(strsplit(args,"="))[1]
        val <- unlist(strsplit(args,"="))[2]

        arg <- gsub(" ", "", arg)
        val <- gsub(" ", "", val)

        vals[[cnt]] <- val
        nms[cnt] <- arg
        cnt <- cnt + 1
    }
    names(vals) <- nms
    return(vals)

}

This function returns the dependent variable of a linear regression

getDepVar <- function(call) {
    call<-parseCall(call)
    formula<-call$formula
    out<-unlist(strsplit(formula,"~")[1])
    return(out[1])
    }

And finally, this function does the magic:

ipredict <-function(model,newdata) {
    Y<-getDepVar(model$call)
    P<-predict(model,newdata=newdata)
    for (i in seq(1,dim(newdata)[1])) {
        if (is.na(newdata[i,Y])) {
            newdata[i,Y]<-predict(model,newdata=newdata[1:i,])[i]
            P[i]<-newdata[i,Y]
        }
    }
    return(P)
}

Example usage (based on my question):

#A function to calculate lags
lagmatrix <- function(x,max.lag){embed(c(rep(NA,max.lag),x),max.lag)}
lag <- function(x,lag) {
    out<-lagmatrix(x,lag+1)[,lag]
    return(out[1:length(out)-1])
}

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

#Fit linear model
model<-lm(y~A+B+C+I(lag(y,1)))
summary(model)

#Forecast linear model
A<-c(A,rnorm(50)) #Assume we know 50 future values of A, B, C
B<-c(B,rnorm(50))
C<-c(C,rnorm(50))
y<-c(y,rep(NA,50))

newdata<-as.data.frame(cbind(y,A,B,C))
ipredict(model,newdata=newdata)
$\endgroup$
  • $\begingroup$ your answer indicates that your question is more about R implementation, than about statistics. So you better ask it at stackoverflow.com. $\endgroup$ – mpiktas Jan 26 '11 at 5:07
  • $\begingroup$ what is the purpose of parseCall function? In function getDepVar you call parseCall and then access formula object from parseCall result, why not accessing it directly from the call? When you have formula object you can get the response name using the all.vars function, try all.vars(formula(y~x)[[2]]). $\endgroup$ – mpiktas Jan 26 '11 at 5:19
  • 1
    $\begingroup$ you probably want lag(y,-1), not lag(y,1). Try cbind(y,lag(y,-1),lag(y,1)) to see the difference. $\endgroup$ – mpiktas Jan 26 '11 at 5:20
  • 1
    $\begingroup$ and final comment. R has a powerful formula interface, use it. Look at the code of lm, and functions model.matrix, model.frame, model.response. Using strsplit and parse on formulas and call objects is a bit of trying to fit a square peg into the round hole. $\endgroup$ – mpiktas Jan 26 '11 at 5:22
  • $\begingroup$ @mpiktas: I posted my second draft, what do you think? $\endgroup$ – Zach Jan 31 '11 at 21:08

Not the answer you're looking for? Browse other questions tagged or ask your own question.