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)

 A: 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)

A: 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.
A: 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)

