require(Matching)


#Run replicates

nrep<-1000
ests.mat.match<-matrix(0,nrow=nrep,ncol=6)

for(irep in 1:nrep){
	#Predictors
	X<-mvrnorm(n,mu=rep(0,p),Sigma=Sig)

	#Data simulation
	Xmat<-cbind(rep(1,n),X)
	eta.Z<-Xmat %*% (c(1,al.select)*al)
	p.Z<-expit(eta.Z)
	Z<-rbinom(n,1,p.Z)
	eta.Y<-Z*theta+Xmat %*% (c(1,be.select)*be)
	
	Y<-eta.Y+rnorm(n)*sig.Y
	#boxplot(Y~Z)

	##########################################################

	#Matching using: all predictors

	p1<-glm(Z~X,family=binomial)
	ps.vals<-fitted(p1)

	rr  <- Match(Y=Y, Tr=Z, X=X, M=1, estimand='ATE');
	ests.mat.match[irep,1]<-rr$est

	rr  <- Match(Y=Y, Tr=Z, X=ps.vals, M=1, estimand='ATE');
	ests.mat.match[irep,2]<-rr$est

	rr  <- Match(Y=Y, Tr=Z, X=ps.vals, M=5, estimand='ATE');
	ests.mat.match[irep,3]<-rr$est
	
	rr  <- Match(Y=Y, Tr=Z, X=ps.vals, M=10, estimand='ATE');
	ests.mat.match[irep,4]<-rr$est

	p1<-glm(Z~X[,confound.list],family=binomial)
	ps.vals<-fitted(p1)

	rr  <- Match(Y=Y, Tr=Z, X=X[,confound.list], M=1, estimand='ATE');
	ests.mat.match[irep,5]<-rr$est

	rr  <- Match(Y=Y, Tr=Z, X=ps.vals, M=1, estimand='ATE');
	ests.mat.match[irep,6]<-rr$est


}

mse.vals<-n*apply((ests.mat.match-theta)^2,2,mean)
name.list<-c('X','M=1','M=5','M=10','Conf.X','Conf.PS')
boxplot(ests.mat.match,names=name.list,ylim=range(-2,6),main='Matching Estimates')
abline(h=theta,col='red')
text(1:ncol(ests.mat.match),rep(6,ncol(ests.mat.match)),format(mse.vals,digits=4))

fname<-paste(filecode,'-PS-Match.pdf',sep='')
dev.print(device = pdf, file=fname,width=11)

boxplot(Y~Z,names=c('Untreated','Treated'))
mean(Y[Z==1])-mean(Y[Z==0])