Rolling analysis with out-of sample I have a model that looks like 
lm(y ~ lag(x, -1) + lag(z, -1))

So basically, this is a time series regression with exogenous variables, and I want to carry out a rolling analysis of sample forecasts, meaning that:
I first used a subsample (e.g., 1990-1995) for estimation, then I performed a one step ahead forecast, then I added one observation and made another one step ahead forecast, and so on.
I have tried to work with rollapply, defining the model as arima(0,0,0) with xreg=lags of the other variables, but that doesn't work. 
Your help would be much appreciated!
 A: Here's a brute-force method, which in general I prefer if a) I can't find an appropriate R function in about 3 minutes, and b) I can see that the brute force function's going to be easy to write.
First, I would start by realigning the variables in a data frame so you don't need to use the lag function:
N <- nrow(y)
df <- as.data.frame(cbind(y[2:N],x[1:(N-1)],z[1:(N-1)]))
colnames(df) <- c("y","x","z")

Then, define a prediction function:
# Start with M observations, gather 1-step-ahead predictions
predict.1 <- function(f, df, M)
{
  P <- nrow(df) - M
  results <- rep(0, P)

  for (i in 1:P) {
    df.pred <- df[M+i,]
    df.est <- df[1:(M+i-1),]
    results[i] <- predict(lm(f, data=df.est), newdata=df.pred)
  }
  results
}

Of course, you could make it more terse, but I'm trying to make it (a little) clearer than it could be.  Using this function looks like:
> # Create sample data
> # Pretend we've "realigned" lagged variables so we don't need to refer to them as lagged.
> x <- rnorm(50)
> z <- rnorm(50)
> y <- x + z + rnorm(50)
> df <- as.data.frame(cbind(y, x, z))
> colnames(df) <- c("y","x","z")
> 
> pred.vals <- predict.1(y~x+z, df, 40)
> pred.vals
 [1] -0.33967757  2.30165856 -0.40084611  0.31978776 -1.75524544
 [6] -0.21552467  0.09107069  0.53836453  0.19864094  2.09003861
> 

It should be pretty obvious how to change the function to accept a parameter for the forecast horizon.  If you want your one-step-ahead forecasts to always use the same number of data points in the history (instead of growing the number of data points in the history by one each step) that's a pretty simple change too.  I pass the formula for the regression function in "f" so that, if I'm comparing different models, I don't have to change the interior of predict.1 each time.
A: I had the same situation and the following R code solved my problem. Although I was using linear regression (lm), you can replace it with ARIMA if you want. 
dd          <- read.csv("SCALE8_data_ready_4_BN_exp1.csv")
windowsSize <- 4000  # training data size
testsize    <- 70    # number of observation to forecast
# load variables from the data set 
delta1   <- dd$delta_t_discrete_3  #$
OSVOS1   <- dd$OSV_OS_discrete     #$
nbAlpha1 <- dd$nb_alpha_corrected  #$
OSVEXT1  <- dd$OSVEXT              #$
h1       <- dd$hidden_state        #$
RMSE     <- matrix(0, 50, 1)
for(k in 0:33)  # run 34 experiments
{
  A         <- k*testsize + 1
  B         <- A + windowsSize - 1
  start_obs <- A
  end_obs   <- B

  delta <- delta1[A:B]
  NbAlpha <- nbAlpha1[A:B]
  OSVOS   <- OSVOS1[A:B]
  OSVEXT  <- OSVEXT1[A:B]

 # ddata  <- data.frame(delta=delta, NbAlpha=NbAlpha, OSVOS=OSVOS, OSVEXT=OSVEXT)
 # output <- paste("Gold_theta0.2per_alpha_0.1_CPTs_from", A, "to", B, 
 #                 "RollingLMTraining.csv")
 # write.csv(ddata, file=output)

  llmm <- lm(OSVEXT~delta + NbAlpha + OSVOS)  # initiate linear regression
  intercept  <- coef(llmm)[1]
  co_delta   <- coef(llmm)[2]
  co_NbAlpha <- coef(llmm)[3]
  co_OSVOS   <- coef(llmm)[4]


  A           <- B + 1
  B           <- B + testsize
  delta       <- delta1[A:B]
  NbAlpha     <- nbAlpha1[A:B]
  OSVOS       <- OSVOS1[A:B]
  OSVEXT      <- OSVEXT1[A:B]
  predict_EXT <- matrix(0, testsize, 1)
  SSE         <- 0
  for(i in 1:testsize)  # do the forecast based on LM results
  {
    predict_EXT[i] <- intercept + delta[i]*co_delta + NbAlpha[i]*co_NbAlpha + 
                      OSVOS[i]*co_OSVOS
    SSE            <- SSE + (predict_EXT[i] - OSVEXT[i])^2
  }
  RMSE[k+1] <- sqrt(SSE/testsize)
 # ddata    <- data.frame(Predicted_EXT=predict_EXT, Real_OSVEXT=OSVEXT)
 # output   <- paste("from", A, "to", B, "RollingLMTesting.csv")
 # write.csv(ddata, file=output)
 print(RMSE[k+1])
}

A: Your question is almost identical to this one. You can use the first step of jbowman answer on this site to make your data.frame and than the rest as I describe in my answer in the linked question with biglm. It is much faster.
