Attempting to use MonteCarlo package in R

Question

I am trying to run a Monte Carlo simulation of a difference in differences estimator, but I am running into an error. Here is the code I am running:

# 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)

And the error that comes up is:

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

This leads me to believe that somewhere something is attempting to access the "0th" element of a vector, which doesn't exist in R as far as I understand. I don't think the part that is doing this arises from my code vs. internal to this package however, and I can't make sense of the code that runs when I run the package.

I am also open to hearing about other methods that will essentially replace simulate() from Stata.

Answer 1

The function passed to MonteCarlo must return a list with named components. Changing line 76 to

return(list("decision" = decision))

should work