以下代码效果很好,但是:我必须在每次运行代码之前更改样本大小 n = 25, 50, ... 和方差估计器。我想用一个循环来解决这个问题。
下面我简单介绍一下代码。在代码中,为给定的样本大小 n 创建了 1000 个回归模型。然后,通过 OLS 估计 1000 个回归模型中的每个回归模型。之后,我根据 1000 个样本中 x3 的不同 beta 值计算 t 统计量。原假设为:H0:beta03 = beta3,即 x3 的计算 beta 值等于我定义为 1 的“真实”值。在最后一步中,我检查了原假设被拒绝的频率(显着性水平 = 0.05)。我的最终目标是创建一个代码,它为每个样本大小和方差估计量吐出原假设的概率拒绝率。如果你们中的任何人可以帮助我,我会很高兴。在这里你可以看到我的代码:
#sample size n = 25, 50, 100, 250, 500, 1000
n <- 50
B <- 1000
#'real' beta values
beta0 <- 1
beta1 <- 1
beta2 <- 1
beta3 <- 1
t.test.values <- rep(NA, B)
#simulation of size
for(rep in 1:B){
#data generation
d1 <- runif(n, 0, 1)
d2 <- rnorm(n, 0, 1)
d3 <- rchisq(n, 1, ncp=0)
x1 <- (1 + d1)
x2 <- (3*d1 + 0.6*d2)
x3 <- (2*d1 + 0.6*d3)
exi <- rchisq(n, 4, ncp = 0)
y <- beta0 + beta1*x1 + beta2*x2 + beta3*x3 + exi
mydata <- data.frame(y, x1, x2, x3)
#ols estimation
lmobj <- lm(y ~ x1 + x2 + x3, mydata)
#extraction
betaestim <- coef(lmobj)[4]
betavar <- vcov(lmobj)[4,4]
#robust variance estimators: hc0, hc1, hc2, hc3
betavar0 <- hccm(lmobj, type="hc0")[4,4]
betavar1 <- hccm(lmobj, type="hc1")[4,4]
betavar2 <- hccm(lmobj, type="hc2")[4,4]
betavar3 <- hccm(lmobj, type="hc3")[4,4]
#t statistic
t.test.values[rep] <- (betaestim - beta3h0)/sqrt(betavar)
}
alpha <- 0.05
test.decision <- abs(t.test.values) < qt(p=c(1-alpha/2), df=n-4)
length(test.decision[test.decision==FALSE])/B