Monte Carlo simulations are useful for studying properties of the sampling distribution of a statistic in a particular controlled environment, where we know the true distribution from which the data are sampled. The statistic of interest could be an estimator, a confidence interval, or a test statistic. Some properties of the statistic that we might be interested in are listed below:
bias of an estimator
variance of an estimator
MSE of an estimator
moments of the distribution of an estimator
quantiles of the distribution of an estimator
coverage rates of a confidence interval
size of a hypothesis test
power function of a hypothesis test
Monte Carlo simulations indicate whether an estimator is consistent, a confidence interval has the correct coverage, or a hypothesis test has the right size (although they do not replace formal proof). We can also use Monte Carlo simulations to compare the biases and MSEs of different estimators or the power curves of various tests.
The idea is that we use computer-generated pseudorandom numbers to create an artificial data set of sample size n to which we apply the statistic of interest. This procedure generates a random draw from the sampling distribution of the statistic.
By repeating the procedure B times independently, we obtain an i.i.d. sample of length B from the sampling distribution of our statistic of interest, which we call a Monte Carlo sample. From the Monte Carlo sample, we can compute empirical counterparts of the features of interest (e.g., bias, MSE, coverage, power, etc.).
8.2 Set up
To set up the Monte Carlo simulation, we need to specify
a sample size n;
a specific parametric distribution F from which we sample our data;
the sampling scheme (i.i.d. or time series process)
the number of Monte Carlo repetitions B.
the Monte Carlo statistic \hat \theta.
For example, if we are interested in the MSE of the sample mean of 100 i.i.d. coin flips, we set n=100, F the Bernoulli distribution with probability 0.5, an i.i.d. sampling scheme, a large number of repetitions (e.g., B=10000), and the sample mean as a Monte Carlo statistic.
If we are interested in the size of a t-test applied to an AR(1) process with parameter 0.8, length 200, and standard normal increments, we set n=200, F the standard normal distribution, the AR(1) process sampling scheme, a large number of replications (e.g., B=10000), and a Monte Carlo statistic \hat \theta containing both the sample mean and the bias-corrected sample variance, since both are needed to compute the t-statistic.
8.3 Monte Carlo algorithm
The Monte Carlo simulation is performed as follows:
Using the specified sampling scheme, draw a sample \{X_1^*, \ldots , X_n^* \} of size n from F using the computer’s random number generator.
Evaluate the statistic \widehat \theta from \{X_1^*, \ldots , X_n^* \}.
Repeat steps 1 and 2 of the experiment B times and collect the estimates from step 2 in the Monte Carlo sample \widehat \theta_{mc} = \{\widehat \theta_1, \ldots, \widehat \theta_B \}.
Evaluate the features of interest from the Monte Carlo sample.
8.4 Evaluate an estimator
If MCsample is the Monte Carlo sample \widehat \theta_{mc} of an estimator \widehat \theta for some parameter \theta (theta), we can evaluate
the MC-estimated sampling mean mean(MCsample):
\widehat \mu_{mc} = \frac{1}{B} \sum_{i=1}^B \widehat \theta_i^*
The MC-estimated p-quantile: quantile(MCsample, p)
We may evaluate whether the estimator is consistent by checking whether the MSE tends to 0 as n gets larger. Let’s try it for the sample mean of an i.i.d. standard normal distributed sample of different sample sizes:
#MC sample meanB=1000#Monte Carlo replicationsgetMCsample=function(n){X=rnorm(n)#standard normal samplemean(X)#sample mean (statistic of interest)}# n=10 (small sample size)MCsample1=sapply(rep(10,B), getMCsample)mse1=var(MCsample1)+(mean(MCsample1)-0)^2mse1
The simulated MSE tends to 0 as n increases. We already derived the theoretical MSE, which is \sigma^2/n. With \sigma^2 = 1, the simulated numbers fit the theoretical ones.
Another example would be checking whether classical standard errors are consistent. That is, se(\overline Y)/sd(\overline Y) \overset{p}{\to} 1 as n \to \infty. We have sd(\overline Y) = \sigma^2/\sqrt n and se(\overline Y) = s_Y/\sqrt n.
#MC standard errorsB=1000#Monte Carlo replicationsgetMCsample=function(n){X=rnorm(n)#standard normal samplese=sd(X)/sqrt(n)#classical standard errorssdt=1/sqrt(n)#true sampling standard deviationse/sdt# return ratio}# n=10 (small sample size)MCsample1=sapply(rep(10,B), getMCsample)mse1=var(MCsample1)+(mean(MCsample1)-1)^2mse1
The MC-estimated MSE of se(\overline Y)/sd(\overline Y) tends to 0, which indicates that the standard error is consistent if the data is sampled from a standard normal distribution.
8.5 Evaluate a confidence interval
An asymptotic 1-\alpha-confidence interval must have a coverage rate that tends to 1-\alpha as the sample size increases.
Let’s check how well a confidence interval with normal quantiles performs when t-quantiles would yield exact confidence intervals. Suppose the data is standard normal with zero mean. The interval has correct coverage if the logical statements lowerbound < 0 and 0 < upperbound are both true. The command (lowerbound < 0) && (0 < upperbound) returns TRUE if both statements are true and FALSE otherwise.
#MC standard errorsB=1000#Monte Carlo replicationsgetMCsample=function(n){X=rnorm(n)#standard normal sample## confidence interval bounds:lowerbound=mean(X)-qnorm(0.975)*sd(X)/sqrt(n)upperbound=mean(X)+qnorm(0.975)*sd(X)/sqrt(n)## check if the true mean 0 is inside the bounds## The command returns TRUE or FALSE(lowerbound<0)&&(0<upperbound)}# n=10 (small sample size)MCsample1=sapply(rep(10,B), getMCsample)sum(MCsample1)/B## relative frequency of TRUEs
A hypothesis test must have asymptotic size \alpha and a power that tends to 1 as the sample size increases.
Let’s evaluate whether the t-test for the mean with \alpha=0.05 works as expected in finite samples if the underlying data has an exponential distribution with parameter \lambda. We consider the two-sided t-test for H_0: \mu = 1. The mean of an exponentially distributed random variable is 1/\lambda, so \mu = 1 if \lambda = 1. If \lambda = 0.8 we have \mu = 1.25.
We consider different sample sizes and evaluate the MC-estimated size for the test under H_0 (i.e., \lambda = 1) and the MC-estimated power of the test for \lambda = 0.8, which is a setting under H_1.
The MC-estimated size is the relative frequency of rejections under H_0, and the MC-estimated power is the relative frequency of rejections under H_1.
#MC standard errorsB=1000#Monte Carlo replicationsgetMCsample=function(n, lambda){X=rexp(n, lambda)# i.i.d. exponential sample## p-value of t-test for mean = 1:pval=t.test(X, mu=1)$p.value## check if p-value is below 0.05 (return TRUE or FALSE):pval<0.05}# n=10 (small sample size)MCsample1.size=sapply(rep(10,B), getMCsample, lambda=1)MCsample1.power=sapply(rep(10,B), getMCsample, lambda=0.8)size=sum(MCsample1.size)/B## relative frequency of TRUEspower=sum(MCsample1.power)/B## relative frequency of TRUEsc(size,power)
[1] 0.099 0.067
# n=100 (moderate sample size)MCsample1.size=sapply(rep(100,B), getMCsample, lambda=1)MCsample1.power=sapply(rep(100,B), getMCsample, lambda=0.8)size=sum(MCsample1.size)/B## relative frequency of TRUEspower=sum(MCsample1.power)/B## relative frequency of TRUEsc(size,power)
[1] 0.052 0.480
# # n=1000 (large sample size)MCsample1.size=sapply(rep(300,B), getMCsample, lambda=1)MCsample1.power=sapply(rep(300,B), getMCsample, lambda=0.8)size=sum(MCsample1.size)/B## relative frequency of TRUEspower=sum(MCsample1.power)/B## relative frequency of TRUEsc(size,power)
[1] 0.050 0.961
The size for n=10 is slightly too large (9.9%) but correct for n=100 and n=300 (5%). For the alternative \mu = 1.25, the test has almost no power for n=10, a power of 48% for n=100, and for n=300, a power of 96%.