Ускорьте случайную цепь Маркова в R с помощью data.table или парелелизации

Я пытаюсь ускорить моделирование методом Монте-Карло дискретной неоднородной по времени цепи Маркова с помощью data.table или некоторой формы распараллеливания. Используя случайные фиктивные матрицы перехода TM, я имитирую временные шаги nSteps в каждом из N симуляций и, начиная с вектора начального состояния initialState, записываю следующее обновленное состояние в currentState. На каждом временном шаге матрица I умножает текущее состояние на матрицу перехода TM.

Код 1 с петлей

nStates <- 5 #number of states
initialState <- c(rep(1/nStates, nStates)) #vector with uniform initial states
nSteps <- 10 #number of time steps
N <- 10000 #number of simulations

ind.arr <- matrix(1:(N*nSteps),ncol=nSteps, byrow=TRUE)
currentState <- vector("list",(N*(nSteps))) #collects the nSteps state vectors for each simulation

system.time(
  for (i in 1:N) {
    TM <- matrix(runif(nStates^2), ncol=nStates) #random transition matrix for each time step and each simulation
    currentState[[(ind.arr[i,1])]] <- initialState %*% (TM / rowSums(TM)) #/rowSums(TM) ensures that TM is a transition matrix
    for (t in 2:nSteps){
      TM <- matrix(runif(nStates^2), ncol=nStates)
      currentState[[(ind.arr[i,t])]] <- currentState[[(ind.arr[i,t-1])]] %*% (TM / rowSums(TM))
    }
  })

Код не очень медленный, но мне интересно, может ли отказ от N-цикла ускорить код. Если я помещу тело N-цикла в функцию

statefun <- function(initialState, nSteps, nStates){
  TM <- matrix(runif(nStates^2), ncol=nStates) #random transition matrix for each time step and each simulation
  currentState <- matrix(rep(NA, nSteps*nStates), ncol=nStates)
  currentState[1,] <- initialState %*% (TM / rowSums(TM)) #/rowSums(TM) ensures that TM is a transition matrix
  for (t in 2:nSteps){
    TM <- matrix(runif(nStates^2), ncol=nStates)
    currentState[t,] <- currentState[t-1,] %*% (TM / rowSums(TM))
  }
  return(currentState)
}

и использую data.table, я получаю ошибку, а не желаемый результат

library(data.table)
system.time(dt <- data.table(i=1:N)[, c("s1", "s2", "s3", "s4", "s5") := list(statefun(initialState, nSteps, nStates)), by=i])

#As each simulation run is independent and the call of statefun is expensive, I was hoping that parallelisation helps to accelerate the code, but trying foreach is actually slower than where I started.  

library(foreach)
system.time(res <- foreach(i=1:N, .combine='c') %do% statefun(initialState, nSteps, nStates))

Я ценю любые комментарии о том, как заставить data.table работать или использовать распараллеливание в этом случае. Большое спасибо, Тим

@ РЕДАКТИРОВАТЬ: этот не принимает десятистрочный вывод вызова функции ...

system.time( #does not work 
  dt <- data.table(i=1:N)[,c("s1", "s2", "s3", "s4", "s5"):=as.list(statefun(initialState, nSteps, nStates)),by=i]
)