# 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 # 1重复n次 X_2 <- rnorm(n, mean = 10, sd = 1.5) # Draw realizations form a normal distribution X_3 <- runif(n, min = 0.2, max = 8) # Draw realizations form a t-distribution X <- cbind(X_1, X_2, X_3) # Save as a Nx3-dimensional data matrix, 合并列；rbind是合并行 beta <- c(1, -5, 5) # c默认列向量 ## Generate realizations from the heteroscadastic error term eps <- 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, 如果不remove，会显示X_1, X_2, X_3 colnames(Var_theo) <- c("", "", "") # remove col name, 如果不remove，会显示X_1, X_2, X_3 ## 如果不remove ## X_1 X_2 X_3 ## X_1 7.22896936 -0.683424680 -0.072931570 ## X_2 -0.68342468 0.069232031 -0.007624144 ## X_3 -0.07293157 -0.007624144 0.057342462 Var_theo ## 7.22896936 -0.683424680 -0.072931570 ## -0.68342468 0.069232031 -0.007624144 ## -0.07293157 -0.007624144 0.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_hat｜x, and Var_beta3_hat｜x ## Solution ## R code for Monte-Carlo Simulation library("sandwich") # HC robust variance estimation MC_reps <- 10000 # Number of Mont-Carlo replications VarHC3_estims <- matrix(NA, 3, MC_reps) # Container to collect the results, matrix(data, nrow, ncol, ...) # 3行因为三种数据，和上面Var_theo行格式一样，10000列用来做模拟 for(r in 1:MC_reps){ ## Generate new relisations from the heteroscedastic error term, 因为conditionally on X，generate new eps eps <- rnorm(n, mean = 0, sd = abs(X_3)) # 10000个eps ## Generate new relizations from the dependent variable，得到new Y Y <- X %*% beta + eps # %*%矩阵乘法 # 10000个Y ## OLS estimation lm_fit <- lm(Y ~ X) ## Robust estimation of the variance of hateta VarHC3_estims[,r] <- diag(sandwich::vcovHC(lm_fit, type = "HC3"))} # 10000个var # 容器是matrix，每次生成的三个元素，以列的形式填充进去 VarHC3_estims_means <- rowMeans(VarHC3_estims) # 返回()内的行均值组成的向量，每种数据的10000个值做平均，得出均值，返回新的向量，返回后默认排列为列向量 ## Compare the theoretical variance Var(hateta_2) and Var(hateta_3) with the means of 10000 variance estimations hat{Var}(hateta_2) and hat{Var}(hateta_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.05945045 plot(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"))