title: How to use Monte Carlo for ARIMA Simulation Function in R
body: <p>Here is the algorithm of what I want to do with R:</p>
  <ol>
  <li>Simulate 10 time series data set from <code>ARIMA</code> model through <code>arima.sim()</code> function</li>
  <li>Split the series into sub-series of possible <code>2s</code>, <code>3s</code>, <code>4s</code>, <code>5s</code>, <code>6s</code>, <code>7s</code>, <code>8s</code>, and <code>9s</code>.</li>
  <li>For each size take a resample the blocks with replacement, for new series and obtain the best <code>ARIMA</code> model from the subseries from each block size through <code>auto.arima()</code> function.</li>
  <li>Obtain for each subseries of each block sizes <code>RMSE</code>.</li>
  </ol>
  <p>The below <code>R</code> function get that done.</p>
  <pre><code>## Load packages and prepare multicore process
  library(forecast)
  library(future.apply)
  plan(multisession)
  library(parallel)
  library(foreach)
  library(doParallel)
  n_cores &lt;- detectCores()
  cl &lt;- makeCluster(n_cores)
  registerDoParallel(cores = detectCores())
  ## simulate ARIMA(1,0, 0)
  #n=10; phi &lt;- 0.6; order &lt;- c(1, 0, 0)
  bootstrap1 &lt;- function(n, phi){
    ts &lt;- arima.sim(n, model = list(ar=phi, order = c(1, 0, 0)), sd = 1)
    ########################################################
    ## create a vector of block sizes
    t &lt;- length(ts)    # the length of the time series
    lb &lt;- seq(n-2)+1   # vector of block sizes to be 1 &lt; l &lt; n (i.e to be between 1 and n exclusively)
    ########################################################
    ## This section create matrix to store block means
    BOOTSTRAP &lt;- matrix(nrow = 1, ncol = length(lb))
    colnames(BOOTSTRAP) &lt;-lb
    ########################################################
    ## This section use foreach function to do detail in the brace
    BOOTSTRAP &lt;- foreach(b = 1:length(lb), .combine = 'cbind') %do%{
      l &lt;- lb[b]# block size at each instance 
      m &lt;- ceiling(t / l)                                 # number of blocks
      blk &lt;- split(ts, rep(1:m, each=l, length.out = t))  # divides the series into blocks
      ######################################################
      res&lt;-sample(blk, replace=T, 10)        # resamples the blocks
      res.unlist &lt;- unlist(res, use.names = FALSE)   # unlist the bootstrap series
      train &lt;- head(res.unlist, round(length(res.unlist) - 10)) # Train set
      test &lt;- tail(res.unlist, length(res.unlist) - length(train)) # Test set
      nfuture &lt;- forecast::forecast(train, model = forecast::auto.arima(train), lambda=0, biasadj=TRUE, h = length(test))$mean        # makes the `forecast of test set
      RMSE &lt;- Metrics::rmse(test, nfuture)      # RETURN RMSE
      BOOTSTRAP[b] &lt;- RMSE
    }
    BOOTSTRAPS &lt;- matrix(BOOTSTRAP, nrow = 1, ncol = length(lb))
    colnames(BOOTSTRAPS) &lt;- lb
    BOOTSTRAPS
    return(list(BOOTSTRAPS))
  }
  </code></pre>
  <p>Calling the function</p>
  <pre><code>bootstrap1(10, 0.6)
  </code></pre>
  <p>I get the below result:</p>
  <pre><code>##              2        3         4        5        6        7         8         9
  ##  [1,] 0.8920703 0.703974 0.6990448 0.714255 1.308236 0.809914 0.5315476 0.8175382
  </code></pre>
  <p>I want to repeat the above <code>step 1</code> to <code>step 4</code> chronologically, then I think of <code>Monte Carlo</code> technology in <code>R</code>. Thus, I load its package and run the below function:</p>
  <pre><code>param_list=list(&quot;n&quot;=10, &quot;phi&quot;=0.6)
  library(MonteCarlo)
  MC_result&lt;-MonteCarlo(func = bootstrap1, nrep=3, param_list = param_list)
  </code></pre>
  <p>expecting to get a like of the below result in <code>matrix</code> form:</p>
  <pre><code>##           [,2]     [,3]      [,4]    [,5]       [,6]      [,7]      [,8]      [,9]
  ##  [1,] 0.8920703 0.703974  0.6990448 0.714255  1.308236  0.809914  0.5315476 0.8175382
  ##  [2,] 0.8909836 0.8457537 1.095148  0.8918468 0.8913282 0.7894167 0.8911484 0.8694729
  ##  [3,] 1.586785  1.224003  1.375026  1.292847  1.437359  1.418744  1.550254  1.30784
  </code></pre>
  <p>but I get the following error message:</p>
  <blockquote>
  <p>Error in MonteCarlo(func = bootstrap1, nrep = 3, param_list = param_list) :
  func has to return a list with named components. Each component has to be scalar.</p>
  </blockquote>
  <p>How can I find my way to obtain a desired result like the above and make the result reproducible?</p>