[SOLVED] CS ## Exercise Chapter 3

## Exercise Chapter 3
## Problem 4

## Background
set.seed(109) # Set the seed of the random number generators
n <- 100 # number of observations## Generate two explanatory variables plus an intercept-variable:X_1 <- rep(1, n) # Intercept # 1nX_2 <- rnorm(n, mean = 10, sd = 1.5) # Draw realizations form a normal distributionX_3 <- runif(n, min = 0.2, max = 8) # Draw realizations form a t-distributionX <- cbind(X_1, X_2, X_3) # Save as a Nx3-dimensional data matrix, rbindbeta <- c(1, -5, 5) # c## Generate realizations from the heteroscadastic error termeps <- rnorm(n, mean = 0, sd = abs(X_3)) # sd = sigma## Dependent variable:Y <- X %*% beta + eps # Y <- X_1*beta_1 + X_2*beta_2 + X_3*beta_3 + eps## Problem 4a compute the theoretical variance of the OLS estimator beta_hat## Solution Var_theo <- solve(t(X) %*% X) %*% t(X) %*% diag(X_3^2) %*% X %*% solve(t(X) %*% X) # Var(beta_hat|X) solve()t()rownames(Var_theo) <- c(“”, “”, “”) # remove row name, removeX_1, X_2, X_3colnames(Var_theo) <- c(“”, “”, “”) # remove col name, removeX_1, X_2, X_3## remove## X_1X_2X_3## X_17.22896936 -0.683424680 -0.072931570## X_2 -0.683424680.069232031 -0.007624144## X_3 -0.07293157 -0.0076241440.057342462Var_theo## 7.22896936 -0.683424680 -0.072931570##-0.683424680.069232031 -0.007624144##-0.07293157 -0.0076241440.057342462## above 2. and 3. diag elements are true variance of beta_hat_2 and beta_hat_3## Problem 4b Monte-Carlo Simulation for Var_beta2_hatx, and Var_beta3_hatx## Solution## R code for Monte-Carlo Simulation library(“sandwich”) # HC robust variance estimationMC_reps <- 10000 # Number of Mont-Carlo replicationsVarHC3_estims <- matrix(NA, 3, MC_reps) # Container to collect the results, matrix(data, nrow, ncol, …) # 3Var_theo10000for(r in 1:MC_reps){## Generate new relisations from the heteroscedastic error term, conditionally on Xgenerate new epseps <- rnorm(n, mean = 0, sd = abs(X_3)) # 10000eps## Generate new relizations from the dependent variablenew YY <- X %*% beta + eps # %*% # 10000Y## OLS estimationlm_fit <- lm(Y ~ X)## Robust estimation of the variance of hatbetaVarHC3_estims[,r] <- diag(sandwich::vcovHC(lm_fit, type = “HC3”))} # 10000var# matrixVarHC3_estims_means <- rowMeans(VarHC3_estims) # ()10000## Compare the theoretical variance Var(hatbeta_2) and Var(hatbeta_3) with the means of 10000 variance estimations hat{Var}(hatbeta_2) and hat{Var}(hatbeta_3)cbind(diag(Var_theo)[c(2,3)], VarHC3_estims_means[c(2,3)])## HC## [,1] [,2]## 0.06923203 0.06578514## 0.05734246 0.05449071## HC3## [,1] [,2]## 0.06923203 0.07258430## 0.05734246 0.05945045plot(x = c(1,2), y = c(0,0), ylim = range(VarHC3_estims[c(2,3),]), type = “n”, axes = FALSE, xlab = “”, ylab = “”)box() # axis(1, c(1,2), labels = c(2,3))axis(2)points(x = rep(1, MC_reps), y = VarHC3_estims[2,], pch = 21, col = gray(.5, .25), bg = gray(.5, .25))points(x = 1, y = VarHC3_estims_means[2], pch = 32, col = “black”, bg = gray(.5, .25))points(x = 1, y = diag(Var_theo)[2], pch = 1)points(x = rep(2, MC_reps), y = VarHC3_estims[3,], pch = 21, col = gray(.5, .25), bg = gray(.5, .25))points(x = 2, y = VarHC3_estims_means[3], pch = 21, col = “black”, bg = “black”)points(x = 2, y = diag(Var_theo)[3], pch = 1)legend(“top”,legend = c(“10000 Variance estimates”, “Mean of variance estimates”, “True Variance”),bty = “n”, pt.bg = c(gray(.5, .75), “black”, “black”), pch = c(21, 21, 1), col = c(gray(.5, .75), “black”, “black”))