尝试在R中使用MonteCarlo包

Question

我试图运行一个差异估计器的蒙特卡洛模拟，但我遇到了一个错误。下面是我正在运行的代码：

# Set the random seed
set.seed(1234567)

library(MonteCarlo)

#Set up problem, doing this before calling the function
# set sample size
n<- 400
# set true parameters: betas and sd of u
b0 <- 1 # intercept for control data (b0 in diffndiff)
b1 <- 1  # shift on both control and treated after treatment (b1 in 
#diffndiff)
b2 <- 2 # difference between intercept on control vs. treated (b2-this is 
#the level difference pre-treatment to compare to coef on treat)
b3 <- 3  # shift after treatment that is only for treated group (b3-this is 
#the coefficient of interest in diffndiff)
b4 <- 0  # parallel time trend (not measured in diffndiff) biases b0,b1 but 
#not b3 that we care about
b5 <- 0 # allows for treated group trend to shift after treatment (0 if 
#parallel trends holds)
su <- 4  # std. dev for errors

dnd <- function(n,b0,b1,b2,b3,b4,b5,su){
#initialize a time vector (set observations equal to n)
timelength = 10
t <- c(1:timelength)


num_obs_per_period = n/timelength #allows for multiple observations in one 
#time period (can simulate multiple states within one group or something)
t0 <- c(1:timelength)
for (p in 1:(num_obs_per_period-1)){
t <- c(t,t0)
}
T<- 5 #set treatment period
g <- t >T
post <- as.numeric(g)

# assign equal amounts of observations to each state to start with (would 
#like to allow selection into treatment at some point)
treat <- vector()
for (m in 1:(round(n/2))){
treat <- c(treat,0)
}
for (m in 1:(round(n/2))){
treat <- c(treat,1)
}
u <- rnorm(n,0,su) #This assumes the mean error is zero

#create my y vector now from the data
y<- b0  + b1*post + b2*treat + b3*treat*post + b4*t + b5*(t-T)*treat*post +u

interaction <- treat*post

#run regression
olsres <- lm(y ~ post + treat + interaction)
olsres$coefficients

# assign the coeeficients
bhat0<- olsres$coefficients[1]
bhat1 <- olsres$coefficients[2]
bhat2<- olsres$coefficients[3]
bhat3<- olsres$coefficients[4]

bhat3_stderr <- coef(summary(olsres))[3, "Std. Error"]

#Here I will use bhat3 to conduct a t-test and determine if this was a pass 
#or a fail
tval <- (bhat3-b3)/ bhat3_stderr

#decision at 5% confidence I believe (False indicates the t-stat was less 
#than 1.96, and we fail to reject the null)
decision <- abs(tval) > 1.96

decision <- unname(decision)
return(list(decision))
} 

#Define a parameter grid to simulate over
from <- -5
to <- 5
increment <- .25
gridparts<- c(from , to , increment)

b5_grid <- seq(from = gridparts[1], to = gridparts[2], by = gridparts[3])

parameter <- list("n" = n, "b0" = b0 , "b1" = b1 ,"b2" = b2 ,"b3" = b3 ,"b4" 
= 
b4 ,"b5" = b5_grid ,"su" = su)

#Now simulate this multiple times in a monte carlo setting
results <- MonteCarlo(func = dnd ,nrep = 100, param_list = parameter)

出现的错误是：

in results[[i]] <- array(NA, dim = c(dim_vec, nrep)) : attempt to select less than one element in integerOneIndex

这让我相信，在某个地方，有人试图访问向量的“0”元素，据我所知，这在R中是不存在的。我不认为做这件事的部分来自我的代码，而不是这个包的内部，而且我不能理解当我运行这个包时运行的代码。

我也乐于听到其他方法，它们将从根本上取代Stata中的simulate()。

Answer 1

传递给MonteCarlo的函数必须返回一个包含命名组件的列表。将第76行更改为

return(list("decision" = decision))

应该行得通