很明显，我之前发布的解决方案是不充分且不优雅的。这是我的第二次尝试，它 100% 解决了我的问题。如果您发现任何错误，请告诉我！如果你们都认为这将是一个更好的地方来获得对我的代码的评论，我将交叉发布到堆栈溢出。

#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)

在其他主题中提出的另一种方法是仅将 arima 函数与 xregs 一起使用。Arima 似乎能够从一组新的 xreg 中做出预测就好了。

好的，我回答了我自己的问题，但我的解决方案可能需要更多测试，并且可能并不完美。建议将不胜感激！

首先，我使用了此处提供的 parseCall 函数的修改版本：

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)

}

此函数返回线性回归的因变量

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

最后，这个函数具有魔力：

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)
}

示例用法（基于我的问题）：

#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)