#enter in the data
y_t <- c(2,6,33,17,2,54)
y_c <- c(1,2,14,14,10,3)

#part a
##Unit level treatment effects
tau_i <- y_t - y_c
##True average treatment effect
tau_bar <- mean(tau_i)

#part b
##True variance
mean_yt <- mean(y_t)
mean_yc <- mean(y_c)
var_yt <- var(y_t)
var_yc <- var(y_c)

var_yt <- sum((y_t - mean_yt)^2) * 1/6
var_yc <- sum((y_c - mean_yc)^2) * 1/6
cov_y <- sum((y_t-mean_yt)*(y_c-mean_yc)) * 1/6


var.usual <- 6/(6-1) * (var_yt/3 + var_yc/3)
var.true <- 6/(6-1) * (var_yt/3 + var_yc/3) + (1/(6-1))*(2*cov_y - var_yt - var_yc)

#Note that the true variance is smaller than the "usual" variance

#part c
exp.fun <- function(y_t, y_c, n.treat){
  #create a treatment vector
   treat.assign <- matrix(0,ncol=1,nrow=length(y_t))
   #Randomly sample and assign treatment
   treat.assign[sample(1:length(treat.assign),n.treat)] <- 1
   #produce observed values
   y_t.obs <- y_t[treat.assign==1]
   y_c.obs <- y_c[treat.assign==0]
   #calculate estimate average treatement effect
   ate <- mean(y_t.obs) - mean(y_c.obs)
   #calculate the standard error
   ate.se <- sqrt(var(y_t.obs)/length(y_t.obs) + var(y_c.obs)/length(y_c.obs))
   #create a list with our desired outputs
   output<- list(y_t.obs = y_t.obs, y_c.obs=y_c.obs, ate = ate, ate.se = ate.se)
   #return the output
   return(output)
}

#Let's make sure the function works
exp.fun(y_t,y_c,3)


#part d
##calculate every possible permutation
perms <- combn(6,3)
#create a matrix to hold all the average treatment effectx
ates <- matrix()
##loop over every combination and calculate the treatment effect
for (i in 1:ncol(perms)){
 ates[i] <-  mean(y_t[perms[,i]]) - mean(y_c[-1* perms[,i]]) 
}
#plot a histogram showing the distribution of the treatment effects
hist(ates,breaks=8, col="darkgreen")