CODEWITHSUNDEEP
rmatrixmontecarlo
diverdan
diverdan

Reputation: 71

How to properly run Monte Carlo simulation using MonteCarlo() in R

I need to run a monte carlo simulation by iterating an experiment a number of times. The functionality of the MonteCarlo package seems to lend itself well to this problem. The experiment involves creating 4 matrices (a, b, c, and d) one after another, for each iteration.

Matrix a:

random.sample <- sample(1:100, 1000, replace = TRUE) 
a = matrix(random.sample, nrow = 10, ncol = 10, byrow = TRUE)

Matrix b: 

b <- apply(a[, 1:10], 2, SMA, n=3)

Matrix c:

c <- pmax(a, b)

Probability sequence:

probability.sequence <- seq(0.1, 0.90, 0.1)

Matrix d:

d <- apply(c, 2, quantile, probs = probability.sequence, na.rm = TRUE)

I would like to iterate the creation of matrices a through d a given amount of times, n. Each iteration should produce a new matrix d containing percentile calculations from the corresponding matrix c. I would like to aggregate the results of the iteration process such that I get a mean, variance, upper confidence limit, and lower confidence limit of the 90% percentile from all matrix d's.

(Approximate) DESIRED OUTPUT:

str(monte.carlo.aggregate.results)
num [1, 1:4] 65 15 85 45 
monte.carlo.aggregate.results         
     mean variance Upper CL Lower CL
[1,]   65       15       85       45

I have read ?MonteCarlo and have struggled to wrap the creation of matrices a through d in a function to be used with the MonteCarlo() function included in the MonteCarlo package. Does anyone have advice for setting up this process for use with MonteCarlo(), or perhaps an alternate solution involving a custom for loop?

Upvotes: 1

Views: 1311

Answers (1)

jay.sf
jay.sf

Reputation: 72593

If you don't insist on using the MonteCarlo package (which is new to me and looks interesting) you could wrap your code into a function.

First, creating a random matrix with only the desired 90%-row. Then, after specifying the amount, you can do the iterations and pack them into a list with lapply(). At the end just create an other matrix with your (rounded) statistics.

Something like this:

set.seed(54897)  # setting seed for sake of reproducibility
probability.sequence <- seq(0.1, 0.90, 0.1)

# iteration function
mx <- function(X, n) {
  random.sample <- sample(1:1e2, 1e3, replace = TRUE) 
  a = matrix(random.sample, nrow = 10, ncol = 10, byrow = TRUE)  
  b <- apply(a[ ,1:10], 2, SMA, n = 3)
  c <- pmax(a, b)
  d <- apply(c, 2, quantile, probs = probability.sequence,  na.rm = TRUE)
  return(d[9, ])  # returns only the desired 90 % row
}

# amount of iterations
amount <- 1e3 

# iteration
mx.lst <- lapply(1:amount, mx)

# matrix for statistics
r = matrix(NA, nrow = 1, ncol = 4, 
           dimnames = list(NULL, c("mean", "variance", "Upper CL", "Lower CL")))
r[, 1] <- (mean(unlist(lapply(mx.lst, mean))))  # mean
r[, 2]  <- sqrt(var(unlist(lapply(mx.lst, mean))))  # variance
r[, 3]  <- r[, 1] - 1.96 * sqrt(r[, 2])  # lower CI 95 %
r[, 4]  <- r[, 1] + 1.96 * sqrt(r[, 2])  # upper CI 95 %

# output
(monte.carlo.aggregate.results <- round(r, 0))
#      mean variance Upper CL Lower CL
# [1,]   82        3       78       86

str(monte.carlo.aggregate.results)
# num 1, 1:4] 82 3 78 86
# - attr(*, "dimnames")=List of 2
# ..$ : NULL
# ..$ : chr [1:4] "mean" "variance" "Upper CL" "Lower CL"

# speed test
amount <- 1e4
system.time(mx.lst <- lapply(1:amount, mx))
#     user      system     elapsed 
#    48.21        0.00       48.31

Note: please check yourself if the used formula for the confidence interval fits to your needs.

Upvotes: 2

Related Questions