# The following code should generate all of the figures
# (with the exception of figure 3) that are included 
# in Handcock & Aldrich (2002).  The full report can be 
# obtained at http://www.csss.washington.edu/Papers

library(reldist)
data(nls)

#
# Now do the RD with smooth=0.4
#
pdf("SmoothFig2.pdf", width=4.0,height=7.5,horiz=FALSE)
par(err=-1)
par(mfrow=c(2,1))
#par(mar=c(2.5,4.1,2.1,2.1), mgp=c(1.5,1,0))
#
fig2a <- reldist(y=recent$chpermwage,yo=original$chpermwage,ci=FALSE,
#        smooth=0.4,
	yowgt=original$wgt,ywgt=recent$wgt,
	cdfplot=TRUE,
	yolabs=seq(-1,3,by=0.5),
	ylabs=seq(-1,3,by=0.5),
	cex=0.8,
        ylab="proportion of the recent cohort",
        xlab="proportion of the original cohort")
title(main="(a)",cex=0.6)
#
fig2b <- reldist(y=recent$chpermwage,yo=original$chpermwage,ci=FALSE,
#        smooth=0.4,
	yowgt=original$wgt,ywgt=recent$wgt,bar=TRUE,
	yolabs=seq(-1,3,by=0.5),
	ylim=c(0,2.5),cex=0.8,
        ylab="Relative Density",
        xlab="Proportion of the Original Cohort")
title(main="(b)",cex=0.6)
#
# Now do the RD with smooth=1.2
#
pdf("SmoothFig3.pdf", width=4.0,height=7.5,horiz=FALSE)
par(err=-1)
par(mfrow=c(2,1))
#
fig3a <- reldist(y=recent$chpermwage,yo=original$chpermwage,ci=FALSE,
#       smooth=0.001,
	yowgt=original$wgt,ywgt=recent$wgt,
	cdfplot=TRUE,
	yolabs=seq(-1,3,by=0.5),
	ylabs=seq(-1,3,by=0.5),
	cex=0.8,
        ylab="proportion of the recent cohort",
        xlab="proportion of the original cohort")
title(main="(a)",cex=0.6)
#
fig3b <- reldist(y=recent$chpermwage,yo=original$chpermwage,ci=FALSE,
        smooth=0.3,
	yowgt=original$wgt,ywgt=recent$wgt,bar=TRUE,
	yolabs=seq(-1,3,by=0.5),
	ylim=c(0,2.5),cex=0.8,
        ylab="Relative Density",
        xlab="Proportion of the Original Cohort")
title(main="(b)",cex=0.6)
#
pdf(file = "SmoothFig4.pdf",width=6.5,height=3.0,horizontal=TRUE)
par(err=-1)
par(mfrow=c(1,3))
#
# Overall r.d. of Y to Yo
# This is a measure of the total differences between Y and Yo.
#
# original$chpermwage is the 
# Predicted change in real log wage from age 16 to 34 for Cohort 1 

g10 <- reldist(y=recent$chpermwage, yo=original$chpermwage, ci=FALSE,
	ywgt=recent$wgt, yowgt=original$wgt,
	yolabs=seq(-1,3,by=0.5),
	ylim=c(0.5,3.0),
	bar=TRUE, quiet=FALSE,
	xlab="proportion of the original cohort")
title(main=paste("(a) entropy = ",format(g10$entropy,digits=3)),cex=0.6)
abline(h=1,lty=2)
#
# calculate the r.d. of Y HAD the 
# only difference from Yo been and ADDITIVE location shift.
# This is a measure of the differences between Y and Yo
# DUE to location differences.
#
# additive location effect
#
# RD of Y to YA
#
g1A <- reldist(y=recent$chpermwage, yo=original$chpermwage,
	ywgt=recent$wgt, yowgt=original$wgt,
	show="effect",
	bar=TRUE, quiet=FALSE,
	ylim=c(0.5,3.0), ylab="",
	ci=FALSE,
	yolabs=seq(-1,3,by=0.5),
              xlab="proportion of the original cohort")
title(main=paste("(b) entropy = ",format(entropy(g1A,g10),digits=3)),cex=0.6)
abline(h=1,lty=2)
#
# calculate the r.d. of Y with the same ADDITIVE location as Yo
# to Yo.
# This is a measure of the differences between Y and Yo
# NOT DUE to the ADDITIVE difference in location.
#
# additive shape effect
#
# RD of YA to Yo
#
gA0 <- reldist(y=recent$chpermwage, yo=original$chpermwage, ci=FALSE,
	ywgt=recent$wgt, yowgt=original$wgt,
	show="residual",
	bar=TRUE, quiet=FALSE,
	ylim=c(0.5,3.0), ylab="",
	yolabs=seq(-1,3,by=0.5),
              xlab="proportion of the original cohort")
title(main=paste("(c) entropy = ",format(gA0$entropy,digits=3)),cex=0.6)
abline(h=1,lty=2)
#
# print out polarizations
#
format(rpy(y=recent$chpermwage,yo=original$chpermwage,
       ywgt=recent$wgt,yowgt=original$wgt,pvalue=TRUE),
       digits=3)
format(rpluy(y=recent$chpermwage,yo=original$chpermwage,
       ywgt=recent$wgt,yowgt=original$wgt,pvalue=TRUE),
       digits=3)
format(rpluy(y=recent$chpermwage,yo=original$chpermwage,
       ywgt=recent$wgt,yowgt=original$wgt,upper=TRUE,pvalue=TRUE),
       digits=3)
#
# endeduc is the final education 
#
e1 <- original$endeduc
e1[e1 < 8] <- 8
e1[e1 > 18] <- 18
e2 <- recent$endeduc
e2[e2 < 8] <- 8
e2[e2 > 18] <- 18
pdf("SmoothFig5.pdf", width=4.5,height=4.5,horiz=FALSE)
par(err=-1)
g10 <- rddist(y=e2, yo=e1, pool=1, ci=FALSE, quiet=FALSE,
        ywgt=recent$wgt,yowgt=original$wgt,
              yolabs=sort(unique(e1)),
              ylab="relative density",
              xlab="proportion of the original cohort")
title(sub=paste("entropy = ",format(entropy(g10),digits=3)))
abline(h=1,lty=2)
#
pdf(file = "SmoothFig6.pdf",width=6.5,height=3.0,horizontal=TRUE)
par(err=-1)
par(mfrow=c(1,3))
#
# Overall r.d. of Y to Yo
# This is a measure of the total differences between Y and Yo.
#
# original$chpermwage is the 
# Predicted change in real log wage from age 16 to 34 for Cohort 1 
#
# generate synthetic Recent data in the same age proportions as
# the Original starting ages
#
# Now the Recent compositionally adjusted for attrition
# and age relative to the Original
#
i3x <- sample(seq(along=original$chpermwage), 
        size = 10*length(original$chpermwage), 
        prob=rdsamp(e2,e1,recent$wgt,original$wgt),
		 replace = TRUE)
schpermwage1 <- original$chpermwage[i3x]
wschpermwage1 <- original$wgt[i3x]
#
g10 <- reldist(y=recent$chpermwage, yo=original$chpermwage, ci=FALSE,
	ywgt=recent$wgt, yowgt=original$wgt,
	yolabs=seq(-1,3,by=0.5),
	ylim=c(0.5,3.0),
	bar=TRUE, quiet=FALSE,
	xlab="proportion of the original cohort")
title(main=paste("(a) entropy = ",format(g10$entropy,digits=3)),cex=0.6)
abline(h=1,lty=2)
#
# calculate the r.d. of Y HAD the 
# only difference from Yo been and ADDITIVE location shift.
# This is a measure of the differences between Y and Yo
# DUE to location differences.
#
# additive location effect
#
# RD of Y to YA
#
#par(xaxs="d", yaxs="d")
g1A <- reldist(y=schpermwage1, yo=original$chpermwage,
	yowgt=original$wgt, ywgt=wschpermwage1,
	bar=TRUE, quiet=FALSE,
	ylim=c(0.5,3.0), ylab="",
	ci=FALSE,
	yolabs=seq(-1,3,by=0.5),
              xlab="proportion of the original cohort")
title(main=paste("(b) entropy = ",format(entropy(g1A,g10),digits=3)),cex=0.6)
abline(h=1,lty=2)
#
# calculate the r.d. of Y with the same ADDITIVE location as Yo
# to Yo.
# This is a measure of the differences between Y and Yo
# NOT DUE to the ADDITIVE difference in location.
#
# additive shape effect
#
# RD of YA to Yo
#
gA0 <- reldist(y=recent$chpermwage, yo=schpermwage1, ci=FALSE,
	ywgt=recent$wgt, yowgt=wschpermwage1,
	bar=TRUE, quiet=FALSE,
	ylim=c(0.5,3.0), ylab="",
	yolabs=seq(-1,3,by=0.5),
              xlab="proportion of the original cohort")
title(main=paste("(c) entropy = ",format(gA0$entropy,digits=3)),cex=0.6)
abline(h=1,lty=2)
#
# print out polarizations
#
format(rpy(y=recent$chpermwage,yo=original$chpermwage,
           ywgt=recent$wgt,yowgt=original$wgt,pvalue=TRUE),
       digits=3)
format(rpluy(y=schpermwage1,yo=original$chpermwage,
             ywgt=wschpermwage1,yowgt=original$wgt,
             pvalue=TRUE),
       digits=3)
format(rpluy(y=recent$chpermwage,yo=schpermwage1,
             ywgt=recent$wgt,yowgt=wschpermwage1,
             upper=TRUE, pvalue=TRUE),
       digits=3)
#
# original$chpermwage is the change in perm wages for the Original
# cohort (with nobs > 2 and non-attrited)
#
# el1 is the final education for the Original
# cohort (with nobs > 2 and non-attrited)
# 1= < HS 2= HS 3=HS+ 4= College+
#
# First generate the samples
#
#
el1 <- e1
el1[el1 < 12] <- 1
el1[el1 == 12] <- 2
el1[el1 > 12 & el1 < 16] <- 3
el1[el1 >= 16] <- 4
#
el2 <- e2
el2[el2 < 12] <- 1
el2[el2 == 12] <- 2
el2[el2 > 12 & el2 < 16] <- 3
el2[el2 >= 16] <- 4
#
i3x <- sample(seq(along=original$chpermwage), 
            size = 10*length(original$chpermwage), 
            prob=rdsamp(el2,el1,recent$wgt,original$wgt),
		 replace = TRUE)
schpermwage1 <- original$chpermwage[i3x]
wschpermwage1 <- original$wgt[i3x]
sel1 <- el1[i3x]
#
pwhso  <- schpermwage1[sel1 <= 2 & !is.na(sel1)]
pwsco  <- schpermwage1[sel1 >  2 & !is.na(sel1)]
wgthso <- wschpermwage1[sel1 <= 2 & !is.na(sel1)]
wgtsco <- wschpermwage1[sel1 >  2 & !is.na(sel1)]
#
pwhsr  <- recent$chpermwage[el2 <= 2 & !is.na(el2)]
pwscr  <- recent$chpermwage[el2 >  2 & !is.na(el2)]
wgthsr <- recent$wgt[el2 <= 2 & !is.na(el2)]
wgtscr <- recent$wgt[el2 >  2 & !is.na(el2)]
#
pdf("SmoothFig7.pdf",width=5.5,height=5.5,horiz=TRUE)
par(err=-1)
par(mfrow=c(2,2))
#
spwhso<-sample(pwhso,size=(100000),prob=wgthso/sum(wgthso),replace=TRUE)
spwsco<-sample(pwsco,size=(100000),prob=wgtsco/sum(wgtsco),replace=TRUE)
spwhsr<-sample(pwhsr,size=(100000),prob=wgthsr/sum(wgthsr),replace=TRUE)
spwscr<-sample(pwscr,size=(100000),prob=wgtscr/sum(wgtscr),replace=TRUE)
#
# Try to overlapping densities
#
nbar <- log(length(spwhso), base = 2) + 1
kwidth <- diff(range(spwhso))/nbar * 0.5
#kwidth <- 1.7*1.4*kwidth
kwidth <- 1.4*kwidth
#
dens1 <- density(spwhso, n = 500, width=2.0*kwidth)
plot(x = (dens1$x), y = dens1$y, type = "l", 
  xlab = "change in log permanent wage", ylab = "density", 
  xlim = c(-1, 3),
  ylim=c(0,1.2))
fig7legend <- list(x=c(0.9,0.9),y=c(1.25,1.25))                  
legend(fig7legend,lty=1:2,cex=0.5, bty="n",
legend=c("original cohort","recent cohort"))       
title(main=paste("(a) high-school or less"),cex=0.6)
#
dens2 <- density(spwhsr, n = 500, width=1.7*kwidth)
lines(x = (dens2$x), y = dens2$y, type = "l",lty=2) 
#
# First the two total r.d. of HS to C between cohort
#
g10hs <- reldist(y=pwhsr, yo=pwhso, ci=FALSE,
	ywgt=wgthsr, yowgt=wgthso,
	bar=TRUE, quiet=FALSE,
	ylim=c(0,4),
        xlab="proportion of the original cohort")
title(main=paste("(b) entropy = ",format(g10hs$entropy,digits=3)),cex=0.6)
abline(h=1,lty=2)
#
nbar <- log(length(spwscr), base = 2) + 1
kwidth <- diff(range(spwscr))/nbar * 0.5
#kwidth <- 1.75*1.2*kwidth
kwidth <- 1.2*kwidth
#
dens1 <- density(spwsco, n = 500, width=1.5*kwidth)
plot(x = (dens1$x), y = dens1$y, type = "l", 
  xlab = "change in log permanent wage", ylab = "density", 
  xlim = c(-1, 3),
  ylim=c(0,1.2))
fig7legend <- list(x=c(0.9,0.9),y=c(1.25,1.25))                  
legend(fig7legend,lty=1:2,cex=0.5, bty="n",
legend=c("original cohort","recent cohort"))       
title(main=paste("(c) more than high school"),cex=0.6)
#
dens2 <- density(spwscr, n = 500, width=2*kwidth)
lines(x = (dens2$x), y = dens2$y, type = "l",lty=2) 
#
g10sc <- reldist(y=pwscr, yo=pwsco, ci=FALSE, 
	ywgt=wgtscr, yowgt=wgtsco,
	bar=TRUE, quiet=FALSE,
	ylim=c(0,4),
        xlab="proportion of the original cohort")
title(main=paste("(d) entropy = ",format(g10sc$entropy,digits=3)),cex=0.6)
abline(h=1,lty=2)
#
format(100*cbind(g10hs$deciles,g10sc$deciles,
                 g10sc$deciles-g10hs$deciles),digits=3)
#
collectresults <- c(
  wtd.median(pwhsr,weight=wgthsr),
  wtd.median(pwhso,weight=wgthso),
  wtd.median(pwhsr,weight=wgthsr)-wtd.median(pwhso,weight=wgthso),
  exp(wtd.median(pwhsr,weight=wgthsr)-wtd.median(pwhso,weight=wgthso)),
  g10hs$entropy,
  g10hs$entropy-reldist(y=pwhsr, yo=pwhso, ywgt=wgthsr, yowgt=wgthso, quiet=FALSE,
                     show="residual",graph=FALSE)$entropy,
  reldist(y=pwhsr, yo=pwhso, ywgt=wgthsr, yowgt=wgthso, quiet=FALSE,
                     show="residual",graph=FALSE)$entropy,
  g10hs$rp[2],
  g10hs$rpl[2],
  g10hs$rpu[2]
)
#
collectresults <- cbind(collectresults,c(
  wtd.median(pwscr,weight=wgtscr),
  wtd.median(pwsco,weight=wgtsco),
  wtd.median(pwscr,weight=wgtscr)-wtd.median(pwsco,weight=wgtsco),
  exp(wtd.median(pwscr,weight=wgtscr)-wtd.median(pwsco,weight=wgtsco)),
  g10sc$entropy,
  g10sc$entropy-reldist(y=pwscr, yo=pwsco, ywgt=wgtscr, yowgt=wgtsco, quiet=FALSE,
                     show="residual",graph=FALSE)$entropy,
  reldist(y=pwscr, yo=pwsco, ywgt=wgtscr, yowgt=wgtsco, quiet=FALSE,
                     show="residual",graph=FALSE)$entropy,
  g10sc$rp[2],
  g10sc$rpl[2],
  g10sc$rpu[2])
)
#
dimnames(collectresults)<-list(
  c("median log recent","median log original",
    "diff median log", "median ratio (recent/original)",
    "entropy", "median entropy","shape entropy","MRP","LRP","URP"),
  c("high-school or less","some college")
)
#
print(collectresults)
binn <- 10
pdf(file = "SmoothFig8.pdf",width=6.5,height=3.5,horizontal=TRUE)
par(err=-1)
par(mfrow=c(1,2))
#
rdhsrscr <- rdeciles(y=pwhsr, yo=pwscr, ywgt=wgthsr, yowgt=wgtscr, binn=binn)
rdhsosco <- rdeciles(y=pwhso, yo=pwsco, ywgt=wgthso, yowgt=wgtsco, binn=binn)
#
mscrdhsrscr <- rdeciles(
 y=pwhsr - wtd.median(pwhsr, weight=wgthsr) + wtd.median(pwhso, weight=wgthso),
yo=pwscr - wtd.median(pwscr, weight=wgtscr) + wtd.median(pwsco, weight=wgtsco),
ywgt=wgthsr, yowgt=wgtscr, binn=binn)
#
mhsrdhsrscr <- rdeciles(
 y=pwhso - wtd.median(pwhso, weight=wgthso) + wtd.median(pwhsr, weight=wgthsr),
yo=pwsco - wtd.median(pwsco, weight=wgtsco) + wtd.median(pwscr, weight=wgtscr),
ywgt=wgthso, yowgt=wgtsco, binn=binn)
#
m1rdhsrscr <- rdeciles(yo=pwsco, 
y=pwhsr - wtd.median(pwhsr, weight=wgthsr) + wtd.median(pwhso, weight=wgthso),
yowgt=wgtsco, ywgt=wgthsr, binn=binn)
#
m2rdhsrscr <- rdeciles(y=pwhso, 
yo=pwscr - wtd.median(pwscr, weight=wgtscr) + wtd.median(pwsco, weight=wgtsco),
ywgt=wgthso, yowgt=wgtscr, binn=binn)
#
m3rdhsrscr <- rdeciles(y=pwhsr, 
yo=pwsco - wtd.median(pwsco, weight=wgtsco) + wtd.median(pwscr, weight=wgtscr),
yowgt=wgtsco, ywgt=wgtscr, binn=binn)
#
achange <- binn*(rdhsrscr$x - rdhsosco$x)
armeff  <- binn*(mhsrdhsrscr$x - rdhsosco$x)
ahseff  <- binn*(m1rdhsrscr$x - rdhsosco$x)
asceff  <- binn*(m2rdhsrscr$x - rdhsosco$x)
ainteff <- achange - armeff - ahseff - asceff
#
bloc <- (1:binn)-0.5
barplot(height=achange,space=0,width=1,axes=FALSE,
     xlab="Decile",ylab="Percentage Point Change",
     ylim=c(-20.0,25))
axis(1,labels=paste(1:binn),at=bloc)
axis(2,labels=TRUE,at=seq(-20.0,25,length=10))
title(main="(a) Marginal effects",cex=0.6)
lines(y=(armeff),x=bloc,lty=1)
lines(y=(asceff),x=bloc,lty=3)
lines(y=(ahseff),x=bloc,lty=2)
abline(h=seq(-20,25,length=10),lty=2)
points(y=(armeff),x=bloc,pch=16,cex=0.7)
points(y=(asceff),x=bloc,pch=3,cex=0.7)
points(y=(ahseff),x=bloc,pch=1,cex=0.7)
#
fig8legend <- list(x=c(4,4),y=c(24.3,24.3))                  
legend(fig8legend,pch=c(16,1,3),lty=c(1:3),cex=0.5, bty="n",
legend=c("Change in relative median",
         "High-school shape effect","College shape effect"))
#
armeff  <- binn*(mhsrdhsrscr$x - rdhsosco$x)
ahseff  <- binn*(m3rdhsrscr$x - mhsrdhsrscr$x)
asceff  <- binn*(rdhsrscr$x - m3rdhsrscr$x)
#
barplot(height=achange,space=0,width=1,axes=FALSE,
     xlab="Decile",ylab="Percentage Point Change",
     ylim=c(-20.0,25))
axis(1,labels=paste(1:binn),at=bloc)
axis(2,labels=TRUE,at=seq(-20.0,25,length=10))
title(main="(b) Sequential effects",cex=0.6)
lines(y=(armeff),x=bloc,lty=1)
lines(y=(asceff),x=bloc,lty=3)
lines(y=(ahseff),x=bloc,lty=2)
abline(h=seq(-20,25,length=10),lty=2)
points(y=(armeff),x=bloc,pch=16,cex=0.7) # filled octagon
points(y=(asceff),x=bloc,pch=3,cex=0.7)  # cross
points(y=(ahseff),x=bloc,pch=1,cex=0.7)  # octagon
#
fig8legend <- list(x=c(4,4),y=c(24.3,24.3))                  
legend(fig8legend,pch=c(16,1,3),lty=c(1:3),cex=0.5, bty="n",
legend=c("Change in relative median",
         "High-school shape effect","College shape effect"))