Hello!

I am working on a program that aims to make randomization inference of four coefficients from a regression. My dataset is a balanced panel that contains information for 13 states and seven years. I have two "treatment groups." The first one is represented by a dummy variable that takes the value of 1 for state i and 0 otherwise and the second variable takes the value of 1 for two other states, different from i. My treatment effects come from the interaction of these two variables with year dummy variables for 2015 and 2016. The four coefficients come from these interactions.

I am estimating the following equation:

yst = state_s + year_t + beta1*treat1*2015_st + beta2*treat1*2016_st+ bet3*treat2*2015_st+beta4*treat2*2016_st + u_st

The purpose of my program is to estimate this equation several times, such that I observe different values for the coefficients beta_1, beta_2, beta_3, and beta_4. The differences in these coefficients will come from randomly assign the states to treatments 1 and 2. However, I am obtaining the results of one single replication, and then the simulation stops after displaying in red font the error .xxxxxx. I would appreciate if someone can explain to me why I cannot obtain more than one set of estimators. Is it because I am using the command "randomselect" to generate the treatment groups?

My code is as follows:

Code:
gen y15 = 0
replace y15 = 1 if year == 2015

gen y16 = 0
replace y16 = 1 if year == 2016


capture prog drop example
program example, rclass 

*Treatment 1
randomselect if state_fips!=0, gen(treat1) select(state) n(1)  /*I am not considering a state for the randomization*/

*Randomly select two bordering states
randomselect if (treat1!=1 & state_fips!=0), gen(treat2) select(state) n(2)
replace treat2 = 0 if treat1 == 1

*Generate treatment variables
gen treat1_15 = treat1*y15
gen treat1_16 = treat1*y16

gen treat2_15 = treat2*y15
gen treat2_16 = treat2*y16

drop treat1 treat2

*Estimation
xtreg y i.year treat1_15 treat1_16 treat2_15 treat2_16 [weight=total], fe vce(cluster state_fips)

capture return scalar t1_15 = _b[treat1_15] 
capture return scalar t1_16 = _b[treat1_16] 
capture return scalar t2_15 = _b[treat2_15] 
capture return scalar t2_16 = _b[treat2_ 16] 

end

simulate t1_15 = r(t1_15) t1_16 = r(t1_16) t2_15 = r(t2_15) t2_16 = r(t2_16), saving("coefficients.dta", replace) reps(100): example
Thank you!