將monte carlo p值排列成不同樣本量和方差估計量的矩陣

Question

以下代碼效果很好（基於my previous question）。但我必須在每次運行代碼之前更改方差估計器（ols,hc0,hc1,hc2,hc3）。我想用循環來解決這個問題。

此後，我簡要介紹一下代碼。在代碼中，創建了每個樣本大小的1000個迴歸模型（n = 25, 50, 100, 250, 500, 1000）。然後，每1000個迴歸模型都由OLS估算。之後，我根據1000個樣本中的不同beta值計算t統計量。零假設爲：H0: beta03 = beta3，即計算出的β值爲等於我定義爲1的「真實」值。在最後一步中，我檢查零假設被拒絕的頻率（顯着性水平= 0.05）。我的最終目標是創建一個代碼，對每個樣本量和方差估計量進行零假設的優先排除率。因此，結果應該是一個矩陣，而現在我得到一個向量作爲結果。如果有人能幫助我，我會很高興。在這裏你可以看到我的代碼：

library(car) 
sample_size = c("n=25"=25, "n=50"=50, "n=100"=100, "n=250"=250, "n=500"=500, "n=1000"=1000) 

B <- 1000 
beta0 <- 1 
beta1 <- 1 
beta2 <- 1 
beta3 <- 1 
alpha <- 0.05 

simulation <- function(n, beta3h0){ 
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) 
# homoskedastic error term: exi <- rchisq(n, 4, ncp = 0) 
exi <- sqrt(x3 + 1.6)*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) 
} 
mean(abs(t.test.values) > qt(p=c(1-alpha/2), df=n-4)) 
} 

sapply(sample_size, simulation, beta3h0 = 1) 

Answer 1

你不需要雙重嵌套循環。只要確保你的循環內有一個矩陣。更新當前simulation具有以下：

## set up a matrix 
## replacing `t.test.values <- rep(NA, B)` 
t.test.values <- matrix(nrow = 5, ncol = B) ## 5 estimators 

## update/fill a column 
## replacing `t.test.values[rep] <- (betaestim - beta3h0)/sqrt(betavar)` 
t.test.values[, rep] <- abs(betaestim - beta3h0)/sqrt(c(betavar, betavar0, betavar1, betavar2, betavar3)) 

## row means 
## replacing `mean(abs(t.test.values) > qt(p=c(1-alpha/2), df=n-4))` 
rowMeans(t.test.values > qt(1-alpha/2, n-4)) 

現在，simulation將返回長度5的向量對於每個樣品的大小，t-統計的p值的蒙特卡洛估計返回所有5個方差估計。然後，當您撥打sapply時，得到矩陣結果：

sapply(sample_size, simulation, beta3h0 = 1) 
#  n=25 n=50 n=100 n=250 n=500 n=1000 
#[1,] 0.132 0.237 0.382 0.696 0.917 0.996 
#[2,] 0.198 0.241 0.315 0.574 0.873 0.994 
#[3,] 0.157 0.220 0.299 0.569 0.871 0.994 
#[4,] 0.119 0.173 0.248 0.545 0.859 0.994 
#[5,] 0.065 0.122 0.197 0.510 0.848 0.993