在 R 中模拟两阶段最小二乘回归

Question

我是 R 的新手，必须为 class 完成以下任务：

这个练习展示了弱工具的有限样本问题两阶段最小二乘法中的 ments。让结构方程为 yi = beta1 + beta2*x_i + eta_i

而简化形式的方程是 xi = gamma1 + gamma2*z_i + e_i:

设样本为i.i.d，其中z_i ~ N (0; 1)，且(eta_i e_i) ~ N((0 0) , (1 0.5 1 0.5))

设置真实参数值为beta1 = 1； beta2 = 2, gamma1 = 0.

对于下面的每个n和gamma2组合，请使用标准的两阶段最小二乘估计器来估计beta2（自己写下估计器。）重复此实验 1000 次，这样我们就有 1000 个估计的 beta2。

n = 100 和 gamma2 = 1。
n = 100 和 gamma2 = 0:1.
n = 100 和 gamma2 = 0:01.
n = 1000 和 gamma2 = 0:01.

在每种情况下绘制这些 beta2 的核密度，以便我们可以看到抽样分布。使用ggplot将4个子图合为一张图

我花了几个小时试图弄清楚如何将其转化为工作模拟。

我最后的解决方案是这样的：

#We will simulate data for this exercise

rm(list = ls( ) )

options(max.print=999999)

#4 different groups of parameters, only n and gamma2 change

  #First Estimation: n=100 and gamma2=1


set.seed(1313)
  #We write a function of the regression model, to later replicate that function
  #and receive 1000 simulations of it.
  b2hat_1_f = function(){
  n = 100 
  #b1 is the constant. We therefore assign it a number of values equal to n.
  b1 = matrix(1, nrow = 100)
  b2 = matrix(2, nrow = 2)
  gamma1 = matrix(0, nrow = 100)
  gamma2 = matrix(1, nrow = 2)

  #Generate Data
  #We have to take the square root of the given variance matrix as rnorm in R uses the standard deviation.
  eta_i = rnorm(n, mean = matrix(0, nrow = 2), sd = sqrt(matrix(c(1, 0.5, 0.5, 1), nrow = 2)))
  e_i = rnorm(n, mean = matrix(0, nrow = 2), sd = sqrt(matrix(c(1, 0.5, 0.5, 1), nrow = 2)))
  Z = cbind(1, rnorm(n, mean = 0, sd = 1))
  X = gamma1 + Z %*% gamma2 + e_i
  X_hat = cbind(1,X)
  Y = b1 + X_hat %*% b2 + eta_i

  #Estimate b2_hat
  #Formula for beta2 estimator is beta_hat_IV = (Z'X)^-1*Z'Y
  #This translates to 
  b2hat = solve((t(Z)%*%X_hat), t(Z)%*%Y)
  }

  b2hat_1 = replicate(1000, b2hat_1_f())

  #Plot the kernel density of our estimator
  plot(density(b2hat_1), main="Kernel Density of b2hat_1 Estimator")

  #Save data of estimator
  save(b2hat_1,file="b2hat_1.RData")`

然后我对其他三个案例进行了类似的设置，最后尝试将它们组合成这样的图表：

#Combine 4 subgraphs using ggplot2
  library(ggplot2)
  library(reshape2)

  load("b2hat_comb.RData")

  b2s = b2hat_comb[, c("b2hat_1", "b2hat_2", "b2hat_3", "b2hat_4")]
  print(head(b2s))

  #Put estimators into data.frame format for use of ggplot2
  x1 = data.frame(b2hat_1)
  x2 = data.frame(b2hat_2)
  x3 = data.frame(b2hat_3)
  x4 = data.frame(b2hat_4)

  ggplot() +
    # b2hat_1
    geom_density(data=x1, aes(x=x1),colour="blue", size=1) +

    # b2hat_2
    geom_density(data=x2, aes(x=x2) ,colour="red", size=1) +

    #b2hat_3
    geom_density(data=x3, aes(x=x3) ,colour="red", size=1) +

    #b2hat_4
    geom_density(data=x4, aes(x=x4) ,colour="red", size=1)

现在我知道模型的设置可能是不正确的，并且知道这部分是由于缺乏计量经济学知识和对 r 的误解。我只是不知道如何继续，目前正在接近 "crisis-mode"。如果你们中有人能抽空看看我做错了什么，我将不胜感激，例如 X_hat 的 cbind 选项对我来说似乎不正确。此外，我不确定我是否使用了正确的 beta2 公式。

我知道这个页面也不应该帮助我完成我的学位，但 ggplot2 部分实际上是一个 r 相关的问题，因为我不知道我做错了什么来产生4 个估计器的组合密度曲线。我尝试了不同的代码，但 none 非常有效。我在这里收到错误消息：

Don't know how to automatically pick scale for object of type data.frame. Defaulting to continuous.
Error: Aesthetics must be either length 1 or the same as the data (2): x

提前非常感谢您，对于我发牢骚的问题深表歉意。如有必要，我很乐意澄清。

Answer 1

好的，我假设我正确处理了数据，但您可以随时纠正我。我更改了函数的输出以生成一个 1x2 矩阵，组合复制的矩阵，并使用该矩阵创建一个包含两个变量的数据框。然后我通过赋值添加了第三个变量，代表数据集的名称 (x1)。然后我将数据框传递给 ggplot() 并让它使用 geom_density 透明地绘制变量 x。它被编写为接受所有四个组 (x1:x4) 并且应该为每个组生成图：

b2hat_1_f = function(){
  n = 100 
  #b1 is the constant. We therefore assign it a number of values equal to n.
  b1 = matrix(1, nrow = 100)
  b2 = matrix(2, nrow = 2)
  gamma1 = matrix(0, nrow = 100)
  gamma2 = matrix(1, nrow = 2)

  #Generate Data
  #We have to take the square root of the given variance matrix as rnorm in R uses the standard deviation.
  eta_i = rnorm(n, mean = matrix(0, nrow = 2), sd = sqrt(matrix(c(1, 0.5, 0.5, 1), nrow = 2)))
  e_i = rnorm(n, mean = matrix(0, nrow = 2), sd = sqrt(matrix(c(1, 0.5, 0.5, 1), nrow = 2)))
  Z = cbind(1, rnorm(n, mean = 0, sd = 1))
  X = gamma1 + Z %*% gamma2 + e_i
  X_hat = cbind(1,X)
  Y = b1 + X_hat %*% b2 + eta_i

  #Estimate b2_hat
  #Formula for beta2 estimator is beta_hat_IV = (Z'X)^-1*Z'Y
  #This translates to 
  b2hat = solve((t(Z)%*%X_hat), t(Z)%*%Y)
  m <- t(b2hat)
}

b2hat_1 = replicate(1000, b2hat_1_f())

m <- b2hat_1
m <- matrix(data = m , ncol=2)
df <- data.frame(x = m[,1], y = m[,2])
df$group <- "x1"

##with the 'group' variable, you don't need separate geoms
library(ggplot2)
ggplot() +
  geom_density(data = df, aes(x = x,  fill = group),alpha=0.5)

希望对您有所帮助。

在 R 中模拟两阶段最小二乘回归

Simulating a 2-Stage Least Squares regression in R

simulation

regression

r

ggplot2