tvec <- c("Treatment", "Time Dummy") par(mfrow=c(2,1)) for(i in 1:2){ plot(density(beta[,i]), main=substitute(paste('Density of ', beta[a]), list(a=tvec[i]))) }
Or, really, if the name of your indexes is the name of the beta columns:
par(mfrow=c(4,2)) for(i in 2:8){ plot(density(beta[,i]), main=substitute(paste('Density of ', beta[a]), list(a=colnames(beta)[i]))) }
source share