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("http://cameron.econ.ucdavis.edu/mmabook/mma12p2mslmsm.asc")
u <- data[, 1]
e <- data[, 2]
y <- data[, 3]
numobs <- length(u)
simreps <- 10000

Simulation

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] 1.187978

and

approx_sterror

[1] 0.01676286