rm(list=ls(all=T)) gc() gc() library("matrixStats") library("locpol") library("KernSmooth") library("gumbel") library("QRM") library("foreign") options(width=80) setwd("/Users/okuiryo/project/drcate/replication_files") Ma <- read.dta("cattaneo2_old.dta") Ma$mbsmoke <- as.integer(Ma$mbsmoke) -1 Ma$fbaby <- as.integer(Ma$fbaby) - 1 Ma$prenatal1 <- as.integer(Ma$prenatal1) -1 Ma$mmarried <- as.integer(Ma$mmarried) -1 ## Use only non-hispanic white mothers Ma <- Ma[which(Ma$mrace ==1),] Ma <- Ma[which(Ma$mhisp ==0),] n <- nrow(Ma) attach(Ma) D <- mbsmoke Y <- bweight X <- cbind(mage, prenatal1, alcohol, fbaby, mage^2, medu, mage*prenatal1, mage*fbaby, mage*medu, mage*alcohol, nprenatal, mage*nprenatal, deadkids, mage*deadkids) k <- ncol(X) x1 <- mage b <- 35 a <- 15 xx1 <- seq(a, b, 0.01) Ytreated<-Y[D==1] Xtreated<-X[D==1,] Yuntreated<-Y[D==0] Xuntreated<-X[D==0,] ### confidence levels alpha<-c(0.01,0.05,0.1) #### Propensity pi_hatT<-fitted(glm(D~X, family=binomial("logit"))) #### Regression miu1_hatF<-cbind(1,X)%*%as.matrix(lm(Ytreated~Xtreated)$coef) miu0_hatF<-cbind(1,X)%*%as.matrix(lm(Yuntreated~Xuntreated)$coef) #### Conditional ATE given all Xs psi1_I<-D*Y/pi_hatT-(D-pi_hatT)*miu1_hatF/pi_hatT psi0_I<-(1-D)*Y/(1-pi_hatT)+(D-pi_hatT)*miu0_hatF/(1-pi_hatT) psi_I<-psi1_I-psi0_I ### unconditional ATE ate <- mean(psi_I) vate <- ate*rep(1,length(xx1)) data<-data.frame(cbind(psi_I,x1)) ### Gaussian Kernel ker<-gaussK rk<-computeRK(ker) lambda<- 1- (1/(4*sqrt(pi)))/rk ### Ruppert,Sheather and Wand (1995) hRSW_I<- dpill(x=x1, y=psi_I) h_I<-hRSW_I*n^(1/5)*n^(-2/7) ghat_I<-locPolSmootherC(x=x1,y=psi_I,xeval=xx1,bw=h_I,deg=1,kernel=ker)$beta0 # standard error fX_hat_I <- numeric(length(xx1)) sigmasq_hat_I <- numeric(length(xx1)) s_hat_I <- numeric(length(xx1)) for(i in 1:length(xx1)) { fX_hat_I[i]<-mean(ker((x1-xx1[i])/h_I))/h_I sigmasq_hat_I[i]<-mean((psi_I-ghat_I[i])^2*ker((x1-xx1[i])/h_I))/fX_hat_I[i]/h_I sigmasq_hat_I[i] <- sigmasq_hat_I[i] *n/(n-3*k) s_hat_I[i]<-sqrt(rk*sigmasq_hat_I[i]/fX_hat_I[i]) } an<-sqrt(2*log(h_I^(-1)*(b-a))+2*log(lambda^(1/2)/(2*pi))) critical_I <- numeric(3) criticalGumbel_I <- numeric(3) critical_N <- c(2.58,1.96,1.64) for (j in 1:3) { critical_I[j]<-sqrt(an^2-2*log(log((1-alpha[j])^(-1/2)))) criticalGumbel_I[j] <- an - log(log((1-alpha[j])^(-1/2)))/(an) } sg_hat_I<-s_hat_I/sqrt(n*h_I) cblower <- matrix(0,ncol = 3, nrow = length(xx1)) cbupper <- matrix(0,ncol = 3, nrow = length(xx1)) cbGumbellower <- matrix(0,ncol = 3, nrow = length(xx1)) cbGumbelupper <- matrix(0,ncol = 3, nrow = length(xx1)) cbNlower <- matrix(0,ncol = 3, nrow = length(xx1)) cbNupper <- matrix(0,ncol = 3, nrow = length(xx1)) for (j in 1:3) { cblower[,j] <- ghat_I - critical_I[j]*sg_hat_I cbupper[,j] <- ghat_I + critical_I[j]*sg_hat_I cbNlower[,j] <- ghat_I - critical_N[j]*sg_hat_I cbNupper[,j] <- ghat_I + critical_N[j]*sg_hat_I cbGumbellower[,j] <- ghat_I - criticalGumbel_I[j]*sg_hat_I cbGumbelupper[,j] <- ghat_I + criticalGumbel_I[j]*sg_hat_I } ylower <- -650 yupper <- 450 pdfname <- c("bw_cb1p.pdf", "bw_cb5p.pdf", "bw_cb10p.pdf" ) for (j in 1:3) { pdf(pdfname[j]) plot(xx1, ghat_I, ylim=c(ylower,yupper), xlab ="mother's age", ylab="", type="l", lwd = 2) par(new = T) plot(xx1, cblower[,j], ylim=c(ylower,yupper), xlab ="", ylab="", type="l", lty =2) par(new = T) plot(xx1, cbupper[,j], ylim=c(ylower,yupper), xlab ="",ylab="", type="l", lty=2) par(new = T) plot(xx1, cbNlower[,j], ylim=c(ylower,yupper), xlab ="", ylab="", type="l", lty =3) par(new = T) plot(xx1, cbNupper[,j], ylim=c(ylower,yupper), xlab ="",ylab="", type="l", lty=3) par(new = T) plot(xx1, cbGumbellower[,j], ylim=c(ylower,yupper),xlab ="", ylab="", type="l", lty=1) par(new = T) plot(xx1, cbGumbelupper[,j], ylim=c(ylower,yupper), xlab ="",ylab="", type="l", lty=1) par(new = T) plot(xx1, vate, ylim=c(ylower,yupper), xlab ="",ylab="", type="l", lty=4) legend("topright", c("CATEF", "our CB", "PW CB","Gumbel CB","ATE"), lty = c(1,2,3,1,4), lwd = c(2,1,1,1,1),cex =1.2) dev.off() }