Search code examples
rmontecarlostandard-error

R Monte-Carlo stopping criteria


Could you please help me understand which of the standard error of the mean (sem) values: sem or semm is more suitable for estimating how close Monte-Carlo simulation to real mean is.

I mean do I have to calculate sem using observations, or semm using mean value after each observation?

#some data
    x <- c(10,9,8,7,6,5,4,3,2,1,10,9,8,7,6,5,4,3,2,1,10,9,8,7,6,5,4,3,2,1,10,9,8,7,6,5,4,3,2,1)
    x_means <- c()
    sem <- c()
    semm <- c()

    for(i in 1:length(x))
    {
      x_means <- c(x_means, mean(x[1:i]))
      sem <- c(sem, sd(x)/sqrt(i))
      semm <- c(semm, sd(x_means)/sqrt(i))
    }

I want to use sem value a stopping criteria for Monte-Carlo simulation, but understand should i calculate sem from sample or means of the sample?


Solution

  • sem in your code uses all of the simulated values, so it is not even feasible for the stopping decision. I guess what you meant for sem is

    sd(x[1:i])/sqrt(i)
    

    If that is so, sem is the right choice.

    For repeated random observations of i.i.d. Y_k's, we're interested in E[Y_k]. Obvious estimator for each increment is X_k = (1/k)(Y_1+...+Y_k), and we want to evaluate precision of X_k for each k. An obvious choice is sample standard deviation of X_k, which is sqrt(1/(k-1)*sum(Y_i-X_k)^2). We can implement this as follows.

    y <- NULL
    precision <- 1
    while (precision > 0.01){
      y <- c(y,rnorm(1)) # your own Monte-Carlo here, for this example, I chose trivial one
      precision <- sd(y)/sqrt(length(y)-1)
      if (is.na(precision)) precision <- 1
    }