Monte Carlo Experiments

This brief tutorial will cover the basics of Monte Carlo Simulations and all the tools that are needed to solve PSet 5. I will use as an example the errors-in-variables model in the context of linear regression. In particular, we want to illustrate with our experiments what we know from studying the theory, as well as some other aspects of programming in R that were left unsaid in the previous tutorial, such as *loops*, *probability distributions*, and *graphics* in R.

# Monte Carlo Experiments

Monte Carlo Experiments are a broad class of computational algorithms that rely on repeated random sampling to obtain numerical results. Their essential idea is using randomness to solve problems that might be deterministic in principle. Monte Carlo methods are mainly used in three distinct problem classes:

- Generating draws from a probability distribtion
- Numerical integration
- Optimization

in this tutorial, we will cover the very basics of the first point, i.e. the possibility of producing a supposedly endless flow of random variables (numbers) for well-known or new distributions.

# Probability Distributions

One convenient use of R is to provide a comprehensive set of statistical tables. Functions are

provided to evaluate the cumulative distribution function \(F_X(x)=\mathbf{P}(X\le x)\), the probability density function of \(F_X(\cdot)\), \(f_X(\cdot)\) and the \(t\)-th quantile function,

\(F_X^{-1}(\tau)=\inf\{x: F_X(x)>\tau\}\), and to simulate from the distribution:

Distribution | R name | Additional Arguments |
---|---|---|

Beta | `beta` |
`shape1` , `shape2` , `ncp` |

Binomial | `binom` |
`size` , `prob` |

Cauchy | `cauchy` |
`location` , `scale` |

Chi-Squared | `chisq` |
`df` , `ncp` |

Exponential | `exp` |
`rate` |

F | `f` |
`df1` , `df2` , `ncp` |

Gamma | `gamma` |
`shape` , `scale` |

Geometric | `geom` |
`prob` |

Hypergeometric | `hyper` |
`m` , `n` , `k` |

Log-normal | `lnorm` |
`meanlog` , `sdlog` |

Logistic | `logis` |
`location` , `scale` |

Negative Binomial | `nbinom` |
`size` , `prob` |

Normal | `norm` |
`mean` , `sd` |

Poisson | `pois` |
`lambda` |

Signed Rank | `signrank` |
`n` |

t | `t` |
`df` , `ncp` |

Uniform | `unif` |
`min` , `max` |

Weibull | `weibull` |
`shape` , `scale` |

Wilcoxon | `wilcox` |
`m` , `n` |

Prefix the name given here by `d`

for the **density**, `p`

for the **CDF**, `q`

for the **quantile** function and `r`

for **simulation** (random deviates). Let's illustrate this for the normal distribution. If you're interested in the density, first argument is

x

```
dnorm(x=1.64)
```

```
[1] 0.1039611
```

q for CDF

```
pnorm(q=0)
```

```
[1] 0.5
```

p for quantile function

```
qnorm(p=0.05)
```

```
[1] -1.644854
```

and n for the radom generator function

```
rnorm(n=5)
```

```
[1] -0.1586431 -0.4372741 1.0921215 -0.1313621 -0.7906500
```

except for `rhyper`

, `rsignrank`

and `rwilcox`

, for which it is nn. In not quite all cases is the non-centrality parameter `ncp`

currently available: see the on-line help for details.

# Loops and Conditional Statements

## if Statement

The conditional construction follows the form `if(exp 1){ exp 2 }else{ exp3}`

. For instance

```
# Generate a random number in [0,1]
r.number <- runif(1)
# Loop
if(r.number>0.3){
print("ok number")
}else {
print("small number")
}
```

```
[1] "small number"
```

There is a vectorized version of the if/else construct, the `ifelse`

function. Using previous example,

```
# Loop
ifelse(r.number>0.3,"ok number","small number")
```

```
[1] "small number"
```

## Repetitive execution: for and while

There is also a `for`

loop construction. For a loop variable say `i`

in a vector expression (most frequently, a sequence), an expression is repeatedly evaluated. For instance,

```
for(i in seq(1:5)){
i^2
}
```

Other looping commands are `repeat`

and `while`

.

# The Model

Let's summarize this model by saying the unobservable variable of interest is \(Z^{*}\), and its observed but mismeasured counterpart is

\[

Z=Z^{*}+\eta

\]

where \(\nu\) denotes the measurement error. Classical Measurement Error assumes that 1) \(Z^{*}\) and \(\eta\) are independent, and 2) \(\mathbb{E}(\eta)=0\).

The textbook example (and conclusion) shows that error in the dependent variable is generally ignored in regression analysis because it simply gets absorbed into the error residual. The only sensible effect is an increase in variability. Hence, inference is more conservative but remains valid nonetheless.

For instance, consider underlying DGP (multivariate normal distribution)

\[

W=

\begin{pmatrix}

Y \

X \

\eta

\end{pmatrix}=\mathcal{N}_3\left(\mathbf{0},

\begin{pmatrix}

\sigma^2_Y & \sigma_{Y,X} & 0 \

\sigma_{X,Y} & \sigma^2_X & 0 \

0 & 0 & \sigma^2_{\eta}

\end{pmatrix}\right)

\]

and you observe

\[

\begin{align*}

Y^{*}&=Y+\eta\

X^{*}&=X+\eta\

\end{align*}

\]

Suppose we're interested in parameter

\[

\beta=\frac{\sigma_{Y,X}}{\sigma_{X}^2}=\frac{\mathbf{E}(YX)}{\mathbf{E}(X^2)}

\]

Draw a random sample of size \(n\). We consider estimators

\[

\begin{align*}

\hat{\beta}_1&=\frac{\sum_{i=1}^nX_iY^{*}_i}{\sum_{i=1}^nX_i^2}\:\:\:

\text{ and }\:\:\:

\hat{\beta}_2=\frac{\sum_{i=1}^nY_iX^{*}_i}{\sum_{i=1}^n (X^{*}_i)^2}

\end{align*}

\]

i.e. in model 1, \(Y\) is mismeasured whereas in model 2 variable \(X\) is mismeasured. It's very easy to show that \(\hat{\beta}_1\) is unbiased and consistent estimator (recall Pset 0). However,

\[

\hat{\beta}_2\overset{P}{\rightarrow}\beta\left(\frac{\sigma^2_X}{\sigma^2_X+\sigma^2_{\eta}}\right)\le\beta

\:\:\:\text{ as }\:\:\:n\to\infty

\]

We can illustrate this with a simple Monte Carlo experiment:

- Draw a sample of size \(1000\) from the underlying multivariate normal distribution with parameters \(\sigma^2_y=1\), \(\sigma_x^2=0.5\), \(\sigma_{\eta}^2=2\), and \(\sigma_{y,x}=0.5\).
- Generate \(Y^{*}\) and \(X^{*}\) as described earlier.
- Calculate \(\hat{\beta}_1\) and \(\hat{\beta}_2\).
- Repeat previous steps \(500\) times.

Then the parameters of the DGP

```
mu <- c(0,0,0)
sigma_y <- 1
sigma_x <- 0.5
sigma_e <- 2
rho=0.5
sigma_yxe <- matrix(c(sigma_y,rho,0,rho,sigma_x,0,0,0,sigma_e),ncol=3)
```

Parameter of interest

```
b <- rho*sqrt(sigma_y/sigma_x)
```

Sample size

```
N <- 1000
```

Number of simulations

```
n.sims <- 500
```

Last but not least, the `seed`

. The seed will allow us to replicate and reproduce results. This is one of the most important practices in Monte Carlo Experiments

```
set.seed(1984)
```

We want to store 6 estimates: \(\beta_0\) and \(\beta_1\) for three models:

- Without measurement error
- Only \(Y\) is measured with error.
- Only \(X\) is measured with error.

The store matrix

```
out <- matrix(NA,nrow=n.sims,ncol=6)
```

at each of the 500 replications, we are generating vector \(Z\)

```
# Call the package
require(mvtnorm)
```

```
Loading required package: mvtnorm
```

```
Z <- rmvnorm(n=N,mean = mu, sigma = sigma_yxe)
```

Contaminated data

```
Ye <- Z[,1]+Z[,3]
Xe <- Z[,2]+Z[,3]
```

Estimate the BLP for the three cases

```
fit1 <- lm(Z[,1]~Z[,2])
fit2 <- lm(Ye~Z[,2])
fit3 <- lm(Z[,1]~Xe)
```

Then the 6 coefficients are given by

```
one.sim <- c(fit1$coefficients,fit2$coefficients,fit3$coefficients)
one.sim
```

```
(Intercept) Z[, 2] (Intercept) Z[, 2] (Intercept)
0.009104504 0.988055248 -0.032965339 0.989833654 0.014665308
Xe
0.215282638
```

If we put this in a loop

```
for(i in 1:n.sims){
# Data Generating Process
Z <- rmvnorm(n=N,mean = mu, sigma = sigma_yxe)
# Contaminated Data
Ye <- Z[,1]+Z[,3]
Xe <- Z[,2]+Z[,3]
# Estimating true relationships
fit1 <- lm(Z[,1]~Z[,2])
# Let's estimate the slope coefficient
# of a linear model ye=bxe+epsilon
fit2 <- lm(Ye~Z[,2])
fit3 <- lm(Z[,1]~Xe)
out[i,] <- c(fit1$coefficients,fit2$coefficients,fit3$coefficients)
}
```

Let's create a data frame with the results

```
coeff1 <- data.frame(b1 = out[,2])
coeff1$type <- 'Error-free'
coeff2 <- data.frame(b1 = out[,4])
coeff2$type <- 'Contaminated'
```

and combine into your new data frame `vegLengths`

```
betaEst <- rbind(coeff1, coeff2)
```

Call package `ggplot2`

```
require(ggplot2)
```

Plot the histograms. First, when only \(Y\) is mismeasured

```
ggplot(betaEst, aes(b1, fill = type)) +
stat_bin(binwidth = 0.001) +
geom_histogram(alpha = 0.5, aes(y = ..density..), position = 'identity')
```

```
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
```

and then when only \(X\) is measured with error

```
coeff3 <- data.frame(b1 = out[,6])
coeff3$type <- 'Contaminated'
betaEst <- rbind(coeff1, coeff3)
ggplot(betaEst, aes(b1, fill = type)) +
stat_bin(binwidth = 0.001) +
geom_histogram(alpha = 0.5, aes(y = ..density..), position = 'identity')
```

```
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
```