Math6630_AS_3_Q6_8.R

# Danny Tshitumbu
# ID: 211964756

# Applied Stats 1: Q. 6. 8.

# Standard Monte Carlo
n <- 1e6
X <- rlnorm(n, 0, 1)
e <- rnorm(n, 0, 1)
Y_X <- exp(9 + 2 * log(X) + e)
EY_X <- mean(Y_X)
EY_X

## [1] 98188.84

# Rao-Blackwellized estimator
n2 <- 1e6
E <- rnorm(n2, 0, 1)
EY_Xe <- rep(0, n2)

for (i in 1:n2) {
  X <- rlnorm(1e2, 0, 1)
  EY_Xe[i] <- mean(exp(9 + 2 * log(X) + E[i]))
}
RB_EX_Y <- mean(EY_Xe)
RB_EX_Y

## [1] 98825.15

# Assessing the Performance of the standard Monte Carlo
system.time({
  n <- 1e6
  X <- rlnorm(n, 0, 1)
  e <- rnorm(n, 0, 1)
  Y_X <- exp(9 + 2 * log(X) + e)
  EY_X <- mean(Y_X)
})

##    user  system elapsed 
##   0.390   0.011   0.406

# Assessing the Performance of the Rao-Blackwell method
system.time({
  n2 <- 1e6
  E <- rnorm(n2, 0, 1)
  EY_Xe <- rep(0, n2)
  
  for (i in 1:n2) {
    X <- rlnorm(1e2, 0, 1)
    EY_Xe[i] <- mean(exp(9 + 2 * log(X) + E[i]))
  }
  RB_EX_Y <- mean(EY_Xe)
})

##    user  system elapsed 
##  47.755   0.713  48.906

# Although the Rao-Blackwell method seems to give result close to the analytical
# result (i.e. exp(11.5) = 98715.77), it seems to be slower than the standard Monte Carlo. 
# We also note that when we increase the number of draws for the standard MC to a 
# relatively large numbers (e.g. bigger than 1e6), the standard MC gives a result
# that is closer to the analytical one (i.e. 98715.77).