require(dplyr); require(tibble)
set.seed(1234)
popsize <- 10e6
eps <- rnorm(popsize)
tr <- rbinom(popsize,1,0.5)
C <- 100
wste <- rep(sample(rep(c(1,-1),times=50),100,rep=FALSE),each=popsize/C)
cl <- rep(c(1:100),each=popsize/C)
y <- wste*tr+eps
#create dataframe
dat <- as_tibble(list(y=y,cluster=cl,treat=tr,te=wste,eps=eps))
rm(tr,eps,te,wste,y)
#create sample#
sampsize <- 100000
df <- dat[sample(nrow(dat),sampsize,rep=FALSE),]
df <- arrange(df,cluster)
#number of observations per cluster#
nobscl <- unlist(count(df,cluster)[,2])
#run the regression and compute the correlations#
res <- lm(df$y~df$treat)
df <- add_column(df,residuals=res$residuals)
df <- add_column(df,prod=df$residuals*df$treat)

#de-mean the objects#
meanres <- rep(0,C)
for(i in 1:C){
  meanres[i] <- mean(df$residuals[df$cluster==i])
}

meanreg <- rep(0,C)
for(i in 1:C){
  meanreg[i] <- mean(df$treat[df$cluster==i])
}

meanprod <- rep(0,C)
for(i in 1:C){
  meanprod[i] <- mean(df$prod[df$cluster==i])
}
df <- add_column(df,demres=df$residuals-rep(meanres,times=nobscl))
df <- add_column(df,demreg=df$treat-rep(meanreg,times=nobscl))
df <- add_column(df,proddem=df$prod-rep(meanprod,times=nobscl))

#finally, compute the within-cluster variance according to the paper
rhoepshat <- (var(df$residuals)-var(df$demres))/var(df$residuals)
rhoreghat <- (var(df$treat)-var(df$demreg))/var(df$treat)
rhoepsregha <- (var(df$prod)-var(df$proddem))/var(df$prod)

round(c(rhoepshat,rhoreghat,rhoepsregha),5)
[1] 0.16686 0.00076 0.24875