r-加速时间序列模拟(用于引导)

Question

我需要在具有非标准依赖性的时间序列上运行引导程序。因此，要做到这一点，我需要创建一个函数，通过逐时间调整来模拟时间序列。

testing<-function(){
  sampleData<-as.zoo(data.frame(index=1:1000,vol=(rnorm(1000))^2,x=NA))
  sampleData[,"x"]<-sampleData[,"vol"]+rnorm(1000) #treat this is completely exognenous and unknown in connection to vol
  sampleData<-cbind(sampleData,mean=rollmean(sampleData[,"vol"],k=3,align="right"))
  sampleData<-cbind(sampleData,vol1=lag(sampleData[,"vol"],k=-1),x1=lag(sampleData[,"x"],k=-1),mean1=lag(sampleData[,"mean"],k=-1))
  #get estimate
  mod<-lm(vol~vol1+x1+mean1,data=sampleData)
  res<-mod$residuals
  for(i in 5:1000){
    #recursively estimate
    sampleData[i,"vol"]<-as.numeric(predict(mod,newdata=data.frame(sampleData[i-1,])))+res[i-3]
    #now must update other paramaters
      #first our rolled average
      sampleData[i,"mean"]<-mean(sampleData[(i-3):i,"vol"])
      #reupdate our lagged variables
      sampleData[i,"vol1"]<-sampleData[i-1,"vol"]
      sampleData[i,"mean1"]<-sampleData[i-1,"mean"]
  }
  lm(vol~vol1+x1+mean1,data=sampleData)
}

当我运行此代码并测量运行时间时，我得到

system.time(testing())
user  system elapsed 
2.711   0.201   2.915 

这对我来说是一个小问题，因为我将集成这些代码来构建一个引导程序。这意味着这里所花费的任何时间对于每一步都要乘以大约100。我已经更新了几千次了。这意味着一次跑步需要数小时（到数天）。

有没有办法加快这个代码的速度？

问候，

Matthew

Answer 1

以下是如何避免predict.lm的开销。还要注意的是，我使用了矩阵而不是动物园对象，这会稍微慢一点。您可以看到这在多大程度上减慢了您的代码。这是你为了方便而付出的代价。

testing.jmu <- function() {
  if(!require(xts)) stop("xts package not installed")
  set.seed(21)  # for reproducibility
  sampleData <- .xts(data.frame(vol=(rnorm(1000))^2,x=NA), 1:1000)
  sampleData$x <- sampleData$vol+rnorm(1000)
  sampleData$mean <- rollmean(sampleData$vol, k=3, align="right")
  sampleData$vol1 <- lag(sampleData$vol,k=1)
  sampleData$x1 <- lag(sampleData$x,k=1)
  sampleData$mean1 <- lag(sampleData$mean,k=1)
  sampleMatrix <- na.omit(cbind(as.matrix(sampleData),constant=1))
  mod.fit <- lm.fit(sampleMatrix[,c("constant","vol1","x1","mean1")],
                    sampleMatrix[,"vol"])
  res.fit <- mod.fit$residuals
  for(i in 5:nrow(sampleMatrix)){
    sampleMatrix[i,"vol"] <-
      sum(sampleMatrix[i-1,c("constant","vol1","x1","mean1")] *
          mod.fit$coefficients)+res.fit[i-3]
    sampleMatrix[i,"mean"] <- mean(sampleMatrix[(i-3):i,"vol"])
    sampleMatrix[i,c("vol1","mean1")] <- sampleMatrix[i-1,c("vol","mean")]
  }
  lm.fit(sampleMatrix[,c("constant","vol1","x1","mean1")], sampleMatrix[,"vol"])
}
system.time(out <- testing.jmu())
#    user  system elapsed 
#    0.05    0.00    0.05 
coef(out)
#    constant        vol1          x1       mean1 
#  1.08787779 -0.06487441  0.03416802 -0.02757601

将set.seed(21)调用添加到您的函数中，您会看到我的函数返回的系数与您的相同。