Econometrics and Free Software by Bruno Rodrigues.
RSS feed for blog post updates.
Follow me on Mastodon, twitter, or check out my Github.
Check out my package that adds logging to R functions, {chronicler}.
Or read my free ebooks, to learn some R and build reproducible analytical pipelines..
You can also watch my youtube channel or find the slides to the talks I've given here.
Buy me a coffee, my kids don't let me sleep.

Method of Simulated Moments with R


This document details section 12.5.6. Unobserved Heterogeneity Example. The original source code giving the results from table 12.3 are available from the authors' site here and written for Stata. This is an attempt to translate the code to R.

Consult the original source code if you want to read the authors' comments. If you want the R source code without all the commentaries, grab it here. This is not guaranteed to work, nor to be correct. It could set your pet on fire and/or eat your first born. Use at your own risk. I may, or may not, expand this example. Corrections, constructive criticism are welcome.

The model is the same as the one described here, so I won't go into details. The moment condition used is \( E[(y_i-\theta-u_i)]=0 \), so we can replace the expectation operator by the empirical mean:

\[ \dfrac{1}{N} \sum_{i=1}^N(y_i - \theta - E[u_i])=0 \]

Supposing that \( E[\overline{u}] \) is unknown, we can instead use the method of simulated moments for \( \theta \) defined by:

\[ \dfrac{1}{N} \sum_{i=1}^N(y_i - \theta - \dfrac{1}{S} \sum_{s=1}^S u_i^s)=0 \]

Import the data

You can consult the original source code to see how the authors simulated the data. To get the same results, and verify that I didn't make mistakes I prefer importing their data directly from their website.

data <- read.table("")
u <- data[, 1]
e <- data[, 2]
y <- data[, 3]
numobs <- length(u)
simreps <- 10000


In the code below, we simulate the equation defined above:

usim <- -log(-log(runif(simreps)))
esim <- rnorm(simreps, 0, 1)

isim <- 0
while (isim < simreps) {

    usim = usim - log(-log(runif(simreps)))
    esim = esim + rnorm(simreps, 0, 1)

    isim = isim + 1


usimbar = usim/simreps
esimbar = esim/simreps

theta = y - usimbar - esimbar

theta_msm <- mean(theta)
approx_sterror <- sd(theta)/sqrt(simreps)

These steps yield the following results:

## Theta MSM= 1.188 Approximate Standard Error= 0.01676