Assignment 02

Instructions

1. Fork this repository to your GitHub account.
2. Write your solutions in R Markdown in a file named solutions.Rmd.
3. When you are ready to submit your assignment, initiate a pull request. Title your pull request "Submission".

To update your fork from the upstream repository:

1. On your fork, e.g. https://github.com/jrnold/Assignment_02 click on "New Pull reqest"
2. Set your fork jrnold/Assignment_02 as the base fork on the left, and UW-POLS503/Assignment_02 as the head fork on the right. In both cases the branch will be master. This means, compare any chanes in the head fork that are not in the base fork. You will see differences between the US-POLS503 repo and your fork. Click on "Create Pull Request", and if there are no issues, "Click Merge" A quick way is to use this link, but change the jrnold to your own username: https://github.com/jrnold/Assignment_02/compare/master...UW-POLS503:master.

Libraries used

library("pols503")
library("mvtnorm")
library("ggplot2")
library("dplyr")
library("broom")

If you do not have the pols503 package installed, you can install it with,

library("devtools")
install_github("UW-POLS503/r-pols503")

OLS Estimator

The purpose of this homework is to provide a guided, hands-on tour through the properties of the least squares estimator, especially under common violations of the Gauss Markov assumptions. We will work through a series of programs which use simulated data --- i.e., data created with known properties --- to investigate how these violations affect the accuracy and precision of least squares estimates of slope parameters. Using repeated study of simulated datasets to explore the properties of statistical models is called Monte Carlo experimentation. Although you will not have to write much R code, you will need to read through the provided programs carefully to understand what is happening.

Monte Carlo experiments always produce the same results as analytic proofs for the specific case considered. Each method has advantages and disadvantages: proofs are more general and elegant, but are not always possible. Monte Carlo experiments are much easier to construct and can always be carried out, but findings from these experiments only apply to the specific scenario under study. Where proofs are available, they are generally preferable to Monte Carlo experiments, but proofs of the properties of more complicated models are sometimes impossible or impractically difficult. This is almost always the case for the properties of models applied to small samples of data. Here, we use Monte Carlo not out of necessity but for pedagogical purposes, as a tool to gain a more intuitive and hands-on understanding of least squares and its properties.

All of the simulations in this assignment will follow the same structure:

1. Define a population

2. Repeat $m$ times:

   1. Draw a sample from the population
   2. Run OLS on that sample
   3. Save statistics, e.g. coefficients, standard errors, $p$-values, from the sample regression.

3. Evaluate the distributions of the sample statistics, or summaries thereof, to determine how well OLS recovers the parameters of the population.

In this section, we will work through the code necessary to run a simulation. However, in the problems, functions written for this problem set will do most of the simulation computation. This section is to help you to understand what those functions are doing, and to provide a mapping from the math to the code.

Sampling Distribution of OLS

Start with example in which the population satisfies all the Gauss-Markov assumptions and we run a correctly specified regression on the samples drawn from that population.

In this example, the population model is $$ \begin{aligned}[t] Y_i &= \beta_0 + \sum_{j = 1}^k \beta_j x_{i,j} + \epsilon_i \ \epsilon_i & \sim N(0, \sigma^2) \end{aligned} $$ For a sample $y$ from that population, the OLS regression which will be run is $$ \begin{aligned}[t] y_i &= \hat\beta_0 + \sum_{j = 1}^k \hat\beta_j x_{i,j} + \hat\varepsilon_i \end{aligned} $$ and the estimate of the population variance is $$ \widehat\sigma^2 &= \frac{\sum \hat\varepsilon_i }{n - k - 1} . $$ In this case, the regression run on the samples has the correct specification, but that will not necessarily be true for other exercises.

Let's review the steps that we need to take for this problem.

1. Draw data from the population
2. Estimate the OLS parameters on that data
3. Return data
4. Repeat 1-3 many times

Of these steps, the meat of the problem is in steps 1--3. Once we have that done, we can repeat it many times with a for loop.

The specified population here is $$ Y_i &= \beta_0 + \sum_{j = 1}^k \beta_j x_{i,j} + \varepsilon_i \ \varepsilon_i & \sim N(0, \sigma^2) $$ What is known (or we need to specify) are the regression parameters, $\beta_0, \dots, \beta_j$, the population variance, $\sigma^2$, and the predictors, $\mat{X}$.

Let's choose to simulate from the population, $$ Y_i &= 0 + 1 x_{i,j} + \varepsilon_i $$ This means that $\vec{\beta} = (0, 1)'$

beta <- c(0, 1)

We also need to choose a sample size, $n$, and the values of $\mat{X}$. While OLS in no way requires $\mat{X}$ to be distributed normal, for these simulations it will be convenient to simulate $\mat{X}$ from a normal distribution, since we can easily choose the mean and variance of the sample. Let's specify that that the sample size is $n = 30$, and the mean and standard deviations are $\bar{\vec{x}} = 0$ and $\bar{\vec{x}} = 1$.

n <- 30
x_mean <- 0
x_sd <- 1

Now let's draw a sample of size 30 from a standard normal distribution, $N(0, 1)$, using the function rnorm,

x_raw <- rnorm(n, mean = 0, sd = 1)

Note that even though the sample was drawn from a normal distribution with the mean and standard deviation that we wanted, due to sampling variation, the sample mean and standard deviation will not be exactly 0 and 1, respectively.

mean(x_raw)

## [1] 0.1035479

sd(x_raw)

## [1] 0.8948433

So, substract the sample mean and divide by the sample standard deviation to get a sample with a mean of 0 and standard deviation of 1,

x <- (x_raw - mean(x_raw)) / sd(x_raw)
mean(x)

## [1] -6.531415e-18

sd(x)

## [1] 1

Note, the R function scale will also do this.

We still need to choose a $\sigma^2$. This is somewhat arbitrary, but I'll chose a $\sigma$ that gives an $R^2 = 0.7$. The function ... in pols503

sigma <- 0.01

Finally, create a data frame with these data in them

sim_linear_normal <- function(.data, beta, sigma) {
  # X gives the number of observations in the data
  n <- nrow(.data)
  # Draw data
  # This creates the X matrix
  X <- model.matrix(~., .data)
  print(X)
  # Create E(y | X)
  yhat <- X %*% beta
  # errors drawn from a normal distribution
  epsilon <- rnorm(n, mean = 0, sd = sigma)
  # actual y's
  y <- yhat + epsilon
  .data$y <- y
  # Estimate model
  mod <- lm(y ~ ., data = .data)
  # Return results
  tidy(mod)
}

This function runs a function multiple times:

run_iterations <- function(.iter, FUN, ...) {
  results <- vector(mode = "list", length = .iter)
  p <- progress_estimated(.iter, min_time = 2)
  for (i in seq_len(.iter)) {
    .data <- FUN(...)
    .data[[".iter"]] <- i
    results[[i]] <- .data
    p$tick()$print()
  }
  bind_rows(results)
}

summarize_params <- function(.data) {
  ret <- .data %>%
    group_by(term) %>%
    summarize(estimate_mean = mean(estimate),
              estimate_sd = sd(estimate),
              std_error_mean = mean(std.error),
              estimate_se = sd(estimate) / sqrt(length(estimate)),
              std_error_se = sd(std.error) / sqrt(length(estimate)),
              iter = length(estimate))
  ret
}

Number of observations in the sample:

n <- 10

For data, we need the number of variables, the mean of the variables, standard deviation of the variables, and their correlation,

k <- 2
mu_X <- rep(0, k)
s_X <- rep(1, k)
R_X <- diag(k)
beta <- c(0, rep(1, k))
X <- as.data.frame(rmvnorm(n, mu_X, sdcor2cov(s_X, R_X)))
names(X) <- paste("x", seq_len(k), sep = "")
sigma <- 0.1