In reply to Isengrim.

The first is correct.
MCMean = MC Estimate of mean = avg(PI Estimates) = average of the estimates over all simulations
MCSEM = MC Estimate of std error of mean = std(PI Estimates) = std dev of the estimates over all simulations

The formula that you write as SEM is (close to) the one-sample estimate of the SEM. If you take a single sample from the simulation and compute std(sample)/sqrt(10,000) then you should get a similar estimate for the precision of the mean. In particular, for the example in this article, you can compute
SEM1 = 4*std(Y) / sqrt(nrow(u));
or
SEM2 = 4*std(u2<Y) / sqrt(nrow(u));
for any of the 1000 samples. The 4 is because the function here is 4*mean.

To compute the CIs for a Monte Carlo estimate, use
MCMean +/- 1.96*MCSEM