# Engel curve example of quantile.inf.np CQTE inference for Goldman and Kaplan (2017), "Nonparametric inference on conditional quantile differences and linear combinations, using L-statistics" # Inspired by Deaton and Paxson (1998) # Questions? Comments? kaplandm@missouri.edu SAVE.FLAG <- TRUE # TRUE=save results to file; FALSE=display in console only. tryCatch(source("quantile_inf_np.R"), error=function(w) { cat("First download the file quantile_inf_np.R currently available from http://faculty.missouri.edu/~kaplandm/code/quantile_inf_np.R and save it in the same directory as this file (or elsewhere in R's path)."); stop(w) } ) tryCatch(library(foreign), error=function(w) { cat("Must install package 'foreign' first."); stop(w) } ) NUMYRS <- 12 #five or ten or twelve if (SAVE.FLAG) { OUTSTUB <- paste0(format(Sys.time(),"%Y_%m_%d"), "_np_qte_ex", sprintf("%02dyrs",NUMYRS)) OUTFILE <- paste0(OUTSTUB,".txt") } else { OUTFILE <- "" } # Input UK survey data raw <- vector('list',12) for (i in 12:(12-NUMYRS+1)) { raw[[i]] <- read.dta(sprintf("20%02d_dvhh_ukanon.dta",i),convert.factors=FALSE) } moneycols <- c('P630tp','P601t','P538t','P539t','P604t') keepcols <- c('A062','Gorx','Year',moneycols) for (i in 12:max(10,12-NUMYRS+1)) { names(raw[[i]])[names(raw[[i]]) %in% keepcols] <- tolower(names(raw[[i]])[names(raw[[i]]) %in% keepcols]) } keepcols <- tolower(keepcols); moneycols <- tolower(moneycols) # dats <- vector('list',12) for (i in 12:(12-NUMYRS+1)) { dats[[i]] <- raw[[i]][,keepcols] } # CPI <- c(94.2,95.4,96.7,98.0,100,102.3,104.7,108.5,110.8,114.5,119.6,123.0) #2001-2012; 2005=100 if (NUMYRS>1) { for (i in 11:(12-NUMYRS+1)) { dats[[i]][,moneycols] <- dats[[i]][,moneycols]*CPI[length(CPI)]/CPI[i] } } # dat <- dats[[12]] if (NUMYRS>1) { for (i in 11:(12-NUMYRS+1)) { dat <- rbind(dat,dats[[i]]) } } # f.1m <- (dat\$a062 %in% c(1,3,5)) f.1w <- (dat\$a062 %in% c(2,4,6)) f.1a <- f.1m | f.1w f.2a <- (dat\$a062 %in% c(7:17)) f.0c <- (dat\$a062 %in% c(1:2,7:8,18,23,26)) f.1c <- (dat\$a062 %in% c(3,4,9:10,19,24)) f.2c <- (dat\$a062 %in% c(5,6,11:12,20)) f.2plusc <- (dat\$a062 %in% c(5,6,11:17,20:22,25,27)) f.1a0c <- (dat\$a062 %in% c(1,2)) f.2a0c <- (dat\$a062 %in% c(7,8)) f.3a0c <- (dat\$a062 %in% c(18)) f.4a0c <- (dat\$a062 %in% c(23)) f.5a0c <- (dat\$a062 %in% c(26)) f.1a1c <- (dat\$a062 %in% c(3,4)) f.2a2c <- (dat\$a062 %in% c(11,12)) f.3a3c <- (dat\$a062 %in% c(21)) f.London <- (dat\$gorx==7) f.SE <- (dat\$gorx==8) n.ppl <- 1*(f.1a0c) +2*(f.2a0c+f.1a1c) +3*(f.3a0c) +4*(f.4a0c+f.2a2c) +5*(f.5a0c) +6*f.3a3c n.K <- 1*f.1a1c +2*f.2a2c +3*f.3a3c n.A <- 1*(f.1a0c+f.1a1c) +2*(f.2a0c+f.2a2c) +3*(f.3a0c+f.3a3c) +4*f.4a0c +5*f.5a0c sum(abs(n.ppl-n.K-n.A)) PCE <- dat\$p630tp / n.ppl filter <- (dat\$p630tp!=0 & n.ppl>0) lnPCE.f <- log(PCE[filter]) # p601t: food at home (EFS). p538t: food at and away from home (ONS). p539t: alcohol (ONS). p604t: housing+utilities w.f <- (dat\$p538t[filter]+dat\$p539t[filter])/dat\$p630tp[filter] # all food + alcohol: mixed (some negative w/ children; most near zero; some positive w/ 2-1 median) w.f <- (dat\$p538t[filter]-dat\$p601t[filter])/dat\$p630tp[filter] # takeout/restaurants only (no alcohol): positive, generally (less so p=0.9) w.f <- (dat[["p601t"]][filter]/dat\$p630tp[filter]) # food at home: negative, generally; ~0 for 2-1 median, 2a2c median. # w.str <- "p538t" # like Deaton and Paxson (1998): all food (excl. alc): ~0 median, some negative 0.9 w.str <- "p604t" # housing: big negative w.str <- "p539t" # alcohol only: generally positive. for (w.str in c("p538t","p604t","p539t")) { w.f <- (dat[[w.str]][filter]/dat\$p630tp[filter]) n.K.f <- n.K[filter] n.A.f <- n.A[filter] n.ppl.f <- n.ppl[filter] itreats <- c(2,5) if (w.str=="p538t") itreats <- c(2,5) tmpn <- 10*n.A.f+n.K.f tmph <- hist(tmpn,breaks=10:51-0.5,plot=FALSE) cat(sprintf("obs(%d)=%d; ",tmph\$mids[tmph\$counts!=0], tmph\$counts[tmph\$counts!=0]), file=OUTFILE,append=TRUE,sep='') cat("\n",file=OUTFILE,append=TRUE) ps <- c(0.5,0.75,0.9) #1:3/4 # w <- ifelse(w<0,0,w) treatments <- rbind(c(4,0),c(2,0),c(3,0),c(4,0), c(2,2),c(3,3)) TREATSTRs <- c("4a0c","2a0c","3a0c","4a0c", "2a2c","3a3c") controls <- rbind(c(1,0),c(1,0),c(2,0),c(3,0), c(1,1),c(2,2)) CTLSTRs <- c("1a0c","1a0c","2a0c","3a0c", "1a1c","2a2c") # treatments <- rbind(c(2,0),c(3,0),c(2,2)) # TREATSTRs <- c("2a0c","3a0c","2a2c") # controls <- rbind(c(1,0),c(2,0),c(1,1)) # CTLSTRs <- c("1a0c","2a0c","1a1c") legposs <- rep.int('bottomright',length(CTLSTRs)) if (w.str=="p604t") legposs <- rep.int('topright',length(CTLSTRs)) for (itreat in itreats) { #1:length(TREATSTRs) TREATSTR <- TREATSTRs[itreat] CTLSTR <- CTLSTRs[itreat] treatment <- (n.A.f==treatments[itreat,1] & n.K.f==treatments[itreat,2]) control <- (n.A.f==controls[itreat,1] & n.K.f==controls[itreat,2]) Y <- w.f[treatment | control] X <- lnPCE.f[treatment | control] Tdummy <- as.integer(treatment[treatment | control]) ALPHA <- 0.10 x0s <- quantile(X,p=1:9/10) CI.pt.los <- CI.pt.his <- matrix(NA,length(x0s),length(ps)) for (pind in 1:length(ps)) { ret <- quantile.inf.np(Y=Y, X=cbind(X,Tdummy), x0s=x0s,p=ps[pind],ALPHA=ALPHA, JOINT.FLAG=F,ONESIDED=0, METHOD.TYPE='qte') CI.pt.los[,pind] <- ret[[1]]\$CI.pointwise.lows CI.pt.his[,pind] <- ret[[1]]\$CI.pointwise.highs } if (SAVE.FLAG) pdf(paste0(OUTSTUB,sprintf("_%s_%s_%s_a%02dpt",w.str,TREATSTR,CTLSTR,100*ALPHA),'.pdf'),width=9,height=8,pointsize=12) par(family='serif',mar=c(5.0,6.0,6.0,2.1)) ylim <- c(-0.10,0.05) if (w.str=='p604t' && itreat==2) ylim <- c(-0.20,0.10) plot(range(x0s),range(c(CI.pt.los,CI.pt.his)),type='n', mgp=c(3.5,1,0),cex.main=2,cex.axis=1.5,cex.lab=2, main=sprintf("%s minus %s, Pointwise %d%%, %s", TREATSTR,CTLSTR,100*(1-ALPHA),w.str), xlab="ln(PCE)", ylab="Budget Share Difference", ylim=ylim, xlim=c(4.3,6.3)) lines(c(0,100),c(0,0),col=1,lty=1,lwd=1,type='l') for (pind in 1:length(ps)) { lines(x0s,CI.pt.los[,pind], lwd=3,type='o',pch=pind,lty=2^(pind-1),col=pind) lines(x0s,CI.pt.his[,pind], lwd=3,type='o',pch=pind,lty=2^(pind-1),col=pind) } legend(legposs[itreat],sprintf("%4.2f-quantile",ps), bty='n', inset=c(0,0), cex=1.8, lty=2^(1:length(ps)-1), pch=1:length(ps), col=1:length(ps),lwd=3) if (SAVE.FLAG) dev.off() } # IQR hiQ <- 0.9; loQ <- 0.5; ALPHA <- 0.10 for (itreat in itreats) { TREATSTR <- TREATSTRs[itreat] CTLSTR <- CTLSTRs[itreat] treatment <- (n.A.f==treatments[itreat,1] & n.K.f==treatments[itreat,2]) control <- (n.A.f==controls[itreat,1] & n.K.f==controls[itreat,2]) Y <- w.f[treatment | control] X <- lnPCE.f[treatment | control] Tdummy <- as.integer(treatment[treatment | control]) x0s <- quantile(X,p=1:9/10) ret <- quantile.inf.np(Y=Y, X=cbind(X,Tdummy), x0s=x0s,p=c(hiQ,loQ),ALPHA=ALPHA, PSI=c(1,-1), JOINT.FLAG=F,ONESIDED=0, METHOD.TYPE='qte') if (SAVE.FLAG) pdf(paste0(OUTSTUB,sprintf("_%s_%s_%s_a%02dpt_IQR",w.str,TREATSTR,CTLSTR,100*ALPHA),'.pdf'),width=9,height=8,pointsize=12) par(family='serif',mar=c(5.0,6.0,6.0,2.1)) plot(range(x0s),range(c(CI.pt.los,CI.pt.his)),type='n', mgp=c(3.5,1,0),cex.main=2,cex.axis=1.5,cex.lab=2, main=sprintf("%s minus %s, (%4.2f-%4.2f)-Interquantile Range, Pointwise %d%%, %s", TREATSTR,CTLSTR,hiQ,loQ,100*(1-ALPHA),w.str), xlab="ln(PCE)", ylab="Interquantile Range Difference", ylim=c(-0.10,0.05), xlim=c(4.3,6.3)) lines(c(0,100),c(0,0),col=1,lty=1,lwd=1,type='l') lines(x0s,ret[[1]]\$CI.pointwise.lows, lwd=3,type='o',pch=1,lty=2,col=1) lines(x0s,ret[[1]]\$CI.pointwise.highs, lwd=3,type='o',pch=1,lty=2,col=1) # legend("top",sprintf("(%4.2f minus %4.2f)-quantile difference",hiQ,loQ), # bty='n', inset=c(0,0), cex=1.8, # lty=2, pch=1, col=1,lwd=3) if (SAVE.FLAG) dev.off() } } #END FOR LOOP over w.str