Ускорение моделирования временных рядов (для начальной загрузки)

Мне нужно запустить загрузчик временного ряда с нестандартной зависимостью. Итак, чтобы сделать это, мне нужно создать функцию, которая имитирует временные ряды, делая корректировку времени по времени.

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 для каждого шага. И я обновляю это несколько тысяч раз. Это означает, что один запуск займет от нескольких часов (до нескольких дней).

Есть ли способ ускорить этот код?

С уважением,

Мэтью


r time-series statistics-bootstrap
person MatthewK    schedule 21.08.2012    source источник
comment
Для большего контекста фактическая функция, которую я использую, берет остатки извне и выводит несколько значений (прогноз, параметры). Остатки передаются через tsboot с непараметрической блочной начальной загрузкой. Затем мне нужно повторять это с течением времени, чтобы увидеть, как меняются параметры (и распределения).   -  person MatthewK    schedule 21.08.2012
comment
Будет ли использование sapply делать это быстрее? Как мне заставить sapply получать значения из строки, над которой он сейчас не работает?   -  person MatthewK    schedule 21.08.2012
comment
sapply не поможет. Вам нужно профилировать свой код, чтобы найти узкое место (см. ?Rprof). Переход с зоопарка на xts экономит немного времени, так как немного времени тратится на поднабор. Вы также можете улучшить производительность, избегая накладных расходов на predict.lm, выполняя умножение вручную.   -  person Joshua Ulrich    schedule 21.08.2012
comment
Я думаю, что основным горлышком бутылки является петля   -  person Luciano Selzer    schedule 21.08.2012
comment
@lselzer: это не так; это predict.lm.   -  person Joshua Ulrich    schedule 21.08.2012
comment
@JoshuaUlrich Вы, вероятно, правы, но предсказание почти 1000 раз тоже не помогает.   -  person Luciano Selzer    schedule 21.08.2012
comment
@lselzer: predict — это универсальный метод, который отправляет метод predict.lm.   -  person Joshua Ulrich    schedule 21.08.2012

Ответы (1)


arrow_upward
7
arrow_downward

Вот как избежать накладных расходов 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) к своей функции, и вы увидите, что моя функция возвращает те же коэффициенты, что и ваша.

person Joshua Ulrich    schedule 21.08.2012
comment
@lselzer: Сколько раз я должен говорить, что проблема не в этом? Почему вы ожидаете, что я удалю его, если это не проблема? Не стесняйтесь давать ответ, который не использует цикл... - person Joshua Ulrich; 21.08.2012