## 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:

- Draw from the DGP
- Compute the estimate
- 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