Home > Programming > Monte Carlo simulations using Stata

## Monte Carlo simulations using Stata

Overview

A Monte Carlo simulation (MCS) of an estimator approximates the sampling distribution of an estimator by simulation methods for a particular data-generating process (DGP) and sample size. I use an MCS to learn how well estimation techniques perform for specific DGPs. In this post, I show how to perform an MCS study of an estimator in Stata and how to interpret the results.

Large-sample theory tells us that the sample average is a good estimator for the mean when the true DGP is a random sample from a $$\chi^2$$ distribution with 1 degree of freedom, denoted by $$\chi^2(1)$$. But a friend of mine claims this estimator will not work well for this DGP because the $$\chi^2(1)$$ distribution will produce outliers. In this post, I use an MCS to see if the large-sample theory works well for this DGP in a sample of 500 observations.

A first pass at an MCS

I begin by showing how to draw a random sample of size 500 from a $$\chi^2(1)$$ distribution and how to estimate the mean and a standard error for the mean.

Example 1: The mean of simulated data

. drop _all
. set obs 500
number of observations (_N) was 0, now 500

. set seed 12345
. generate y = rchi2(1)
. mean y

Mean estimation                   Number of obs   =        500

--------------------------------------------------------------
|       Mean   Std. Err.     [95% Conf. Interval]
-------------+------------------------------------------------
y |   .9107644   .0548647      .8029702    1.018559
--------------------------------------------------------------


I specified set seed 12345 to set the seed of the random-number generator so that the results will be reproducible. The sample average estimate of the mean from this random sample is $$0.91$$, and the estimated standard error is $$0.055$$.

If I had many estimates, each from an independently drawn random sample, I could estimate the mean and the standard deviation of the sampling distribution of the estimator. To obtain many estimates, I need to repeat the following process many times:

1. Draw from the DGP
2. Compute the estimate
3. Store the estimate.

I need to know how to store the many estimates to proceed with this process. I also need to know how to repeat the process many times and how to access Stata estimates, but I put these details into appendices I and II, respectively, because many readers are already familiar with these topics and I want to focus on how to store the results from many draws.

I want to put the many estimates someplace where they will become part of a dataset that I can subsequently analyze. I use the commands postfile, post, and postclose to store the estimates in memory and write all the stored estimates out to a dataset when I am done. Example 2 illustrates the process, when there are three draws.

Example 2: Estimated means of three draws

. set seed 12345

. postfile buffer mhat using mcs, replace

. forvalues i=1/3 {
2.         quietly drop _all
3.         quietly set obs 500
4.         quietly generate y = rchi2(1)
5.         quietly mean y
6.         post buffer (_b[y])
7. }

. postclose buffer

. use mcs, clear

. list

+----------+
|     mhat |
|----------|
1. | .9107645 |
2. |  1.03821 |
3. | 1.039254 |
+----------+


The command

postfile buffer mhat using mcs, replace


creates a place in memory called buffer in which I can store the results that will eventually be written out to a dataset. mhat is the name of the variable that will hold the estimates in the new dataset called mcs.dta. The keyword using separates the new variable name from the name of the new dataset. I specified the option replace to replace any previous versions of msc.dta with the one created here.

I used

forvalues i=1/3 {


to repeat the process three times. (See appendix I if you want a refresher on this syntax.) The commands

quietly drop _all
quietly set obs 500
quietly generate y = rchi2(1)
quietly mean y


drop the previous data, draw a sample of size 500 from a $$\chi^2(1)$$ distribution, and estimate the mean. (The quietly before each command suppresses the output.) The command

post buffer (_b[y])


stores the estimated mean for the current draw in buffer for what will be the next observation on mhat. The command

postclose buffer


writes the stuff stored in buffer to the file mcs.dta. The commands

use mcs, clear
list


drop the last $$\chi^2(1)$$ sample from memory, read in the msc dataset, and list out the dataset.

Example 3 below is a modified version of example 2; I increased the number of draws and summarized the results.

Example 3: The mean of 2,000 estimated means

. set seed 12345

. postfile buffer mhat using mcs, replace

. forvalues i=1/2000 {
2.         quietly drop _all
3.         quietly set obs 500
4.         quietly generate y = rchi2(1)
5.         quietly mean y
6.         post buffer (_b[y])
7. }

. postclose buffer

. use mcs, clear

. summarize

Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
mhat |      2,000     1.00017    .0625367   .7792076    1.22256


The average of the $$2,000$$ estimates is an estimator for the mean of the sampling distribution of the estimator, and it is close to the true value of $$1.0$$. The sample standard deviation of the $$2,000$$ estimates is an estimator for the standard deviation of the sampling distribution of the estimator, and it is close to the true value of $$\sqrt{\sigma^2/N}=\sqrt{2/500}\approx 0.0632$$, where $$\sigma^2$$ is the variance of the $$\chi^2(1)$$ random variable.

Including standard errors

The standard error of the estimator reported by mean is an estimate of the standard deviation of the sampling distribution of the estimator. If the large-sample distribution is doing a good job of approximating the sampling distribution of the estimator, the mean of the estimated standard
errors should be close to the sample standard deviation of the many mean estimates.

To compare the standard deviation of the estimates with the mean of the estimated standard errors, I modify example 3 to also store the standard errors.

Example 4: The mean of 2,000 standard errors

. set seed 12345

. postfile buffer mhat sehat using mcs, replace

. forvalues i=1/2000 {
2.         quietly drop _all
3.         quietly set obs 500
4.         quietly generate y = rchi2(1)
5.         quietly mean y
6.         post buffer (_b[y]) (_se[y])
7. }

. postclose buffer

. use mcs, clear

. summarize

Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
mhat |      2,000     1.00017    .0625367   .7792076    1.22256
sehat |      2,000    .0629644    .0051703   .0464698   .0819693


Mechanically, the command

postfile buffer mhat sehat using mcs, replace


makes room in buffer for the new variables mhat and sehat, and

post buffer (_b[y]) (_se[y])


stores each estimated mean in the memory for mhat and each estimated standard error in the memory for sehat. (As in example 3, the command postclose buffer writes what is stored in memory to the new dataset.)

The sample standard deviation of the $$2,000$$ estimates is $$0.0625$$, and it is close to the mean of the $$2,000$$ estimated standard errors, which is $$0.0630$$.

You may be thinking I should have written “very close”, but how close is $$0.0625$$ to $$0.0630$$? Honestly, I cannot tell if these two numbers are sufficiently close to each other because the distance between them does not automatically tell me how reliable the resulting inference will be.

Estimating a rejection rate

In frequentist statistics, we reject a null hypothesis if the p-value is below a specified size. If the large-sample distribution approximates the finite-sample distribution well, the rejection rate of the test against the true null hypothesis should be close to the specified size.

To compare the rejection rate with the size of 5%, I modify example 4 to compute and store an indicator for whether I reject a Wald test against the true null hypothesis. (See appendix III for a discussion of the mechanics.)

Example 5: Estimating the rejection rate

. set seed 12345

. postfile buffer mhat sehat reject using mcs, replace

. forvalues i=1/2000 {
2.         quietly drop _all
3.         quietly set obs 500
4.         quietly generate y = rchi2(1)
5.         quietly mean y
6.         quietly test _b[y]=1
7.         local r = (r(p)<.05)
8.         post buffer (_b[y]) (_se[y]) (r')
9. }

. postclose buffer

. use mcs, clear

. summarize

Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
mhat |      2,000     1.00017    .0625367   .7792076    1.22256
sehat |      2,000    .0629644    .0051703   .0464698   .0819693
reject |      2,000       .0475     .212759          0          1


The rejection rate of $$0.048$$ is very close to the size of $$0.05$$.

Done and undone

In this post, I have shown how to perform an MCS of an estimator in Stata. I discussed the mechanics of using the post commands to store the many estimates and how to interpret the mean of the many estimates and the mean of the many estimated standard errors. I also recommended using an estimated rejection rate to evaluate the usefulness of the large-sample approximation to the sampling distribution of an estimator for a given DGP and sample size.

The example illustrates that the sample average performs as predicted by large-sample theory as an estimator for the mean. This conclusion does not mean that my friend's concerns about outliers were entirely misplaced. Other estimators that are more robust to outliers may have better properties. I plan to illustrate some of the trade-offs in future posts.

Appendix I: Repeating a process many times

This appendix provides a quick introduction to local macros and how to use them to repeat some commands many times; see [P] macro and [P] forvalues for more details.

I can store and access string information in local macros. Below, I store the string hello" in the local macro named value.

local value "hello"


To access the stored information, I adorn the name of the local macro. Specifically, I precede it with the single left quote () and follow it with the single right quote ('). Below, I access and display the value stored in the local macro value.

. display "value'"
hello


I can also store numbers as strings, as follows

. local value "2.134"
. display "value'"
2.134


To repeat some commands many times, I put them in a {\tt forvalues} loop. For example, the code below repeats the display command three times.

. forvalues i=1/3 {
2.    display "i is now i'"
3. }
i is now 1
i is now 2
i is now 3


The above example illustrates that forvalues defines a local macro that takes on each value in the specified list of values. In the above example, the name of the local macro is i, and the specified values are 1/3=$$\{1, 2, 3\}$$.

Appendix II: Accessing estimates

After a Stata estimation command, you can access the point estimate of a parameter named y by typing _b[y], and you can access the estimated standard error by typing _se[y]. The example below illustrates this process.

Example 6: Accessing estimated values

. drop _all

. set obs 500
number of observations (_N) was 0, now 500

. set seed 12345

. generate y = rchi2(1)

. mean y

Mean estimation                   Number of obs   =        500

--------------------------------------------------------------
|       Mean   Std. Err.     [95% Conf. Interval]
-------------+------------------------------------------------
y |   .9107644   .0548647      .8029702    1.018559
--------------------------------------------------------------

. display  _b[y]
.91076444

. display _se[y]
.05486467


Appendix III: Getting a p-value computed by test

This appendix explains the mechanics of creating an indicator for whether a Wald test rejects the null hypothesis at a specific size.

I begin by generating some data and performing a Wald test against the true null hypothesis.

Example 7: Wald test results

. drop _all

. set obs 500
number of observations (_N) was 0, now 500

. set seed 12345

. generate y = rchi2(1)

. mean y

Mean estimation                   Number of obs   =        500

--------------------------------------------------------------
|       Mean   Std. Err.     [95% Conf. Interval]
-------------+------------------------------------------------
y |   .9107644   .0548647      .8029702    1.018559
--------------------------------------------------------------

. test _b[y]=1

( 1)  y = 1

F(  1,   499) =    2.65
Prob > F =    0.1045


The results reported by test are stored in r(). Below, I use return list to see them, type help return list for details.

Example 8: Results stored by test

. return list

scalars:
r(drop) =  0
r(df_r) =  499
r(F) =  2.645393485924886
r(df) =  1
r(p) =  .1044817353734439


The p-value reported by test is stored in r(p). Below, I store a 0/1 indicator for whether the p-value is less than $$0.05|0 in the local macro r. (See appendix II for an introduction to local macros.) I complete the illustration by displaying that the local macro contains the value \(0$$.

. local r = (r(p)<.05)
. display "r'"
0

Categories: Programming Tags: