Calculate expected value of variance using monte carlo simulation

Question

So I have this probability distribution

X = {0      probability 7/8}
      {1/60 probability 1/8}

James his car breaks down N times a year where N ~ Pois(2) and X the repair cost and Y is the total cost caused by James in a year.

I want to calculate the E[Y] and V(Y), which should give me E[X]=15 and V(Y) = 1800

I have this monte Carlo simulation:

expon_dis <- rexp(200, 1/60)

result_matrix2 <- rep(0, 200)
expected_matrix <- rep(0, runs)

for (u in 1:runs){
  expon_dis <- rexp(200, 1/60)
  N <- rpois(200, 2)
  for (l in 1:200){
    result_matrix2[l] <- (expon_dis[l] * (1/8)) * (N[l])
  }
  expected_matrix[u] <- mean(result_matrix2)
}

This code gives the expected value of 15 but the variance is not correct. So what is wrong with this simulation?

Answer 1

Not enough time to read through your code, but i think the error comes with the multiplication.

Below is a very rough implementation, where first you write a function to simulate the cost, given x number of breakdowns:

sim_cost = function(x){
cost = rexp(x,1/60)
prob = sample(c(0,1/60),x,prob=c(7/8,1/8),replace=TRUE)
sum(cost[prob>0])
}

Then generate the number of breakdowns per year:

set.seed(111)
N <- rpois(500000, 2)

Iterate over the years, if 0, we return 0:

set.seed(111)
sim = sapply(N,function(i)if(i==0){0}else{sum(sim_cost(i))})

mean(sim)
[1] 14.98248
var(sim)
[1] 1797.549

You need quite a number of simulations, but above should be a code that you can start to optimize to get it closer.