# Predicting from a simple linear model with lags in R

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)


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)  • @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. – mpiktas Feb 1 '11 at 8:22 • @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? – Zach Feb 1 '11 at 14:21 • @Zach, the problem is with one-step forecast. Apparently it does not work. – mpiktas Feb 1 '11 at 14:30 • @It seems to work with the dyn function, rather than the dynlm function. – Zach Feb 1 '11 at 16:17 • @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. – mpiktas Feb 2 '11 at 7:52

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.

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)