Running 'FigSM.R' [161s]
Running the tests in 'tests/FigSM.R' failed.
Complete output:
  > # 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 SSRN: http://dx.doi.org/10.2139/ssrn.1515775
  > 
  > library(reldist)
  reldist: Relative Distribution Methods
  Version 1.7-2 created on 2023-02-16.
  copyright (c) 2003, Mark S. Handcock, University of California-Los Angeles
   For citation information, type citation("reldist").
   Type help(package="reldist") to get started.
  
  > data(nls, package="reldist")
  > 
  > pdf("Fig1.pdf", width=4.0,height=7.5)
  > par(mfrow=c(2,1))
  > par(err=-1)
  > #par(mar=c(2.5,4.1,2.1,2.1), mgp=c(1.5,1,0))
  > #
  > # Set the random seed
  > #
  > # First original from Splus
  > #
  > #.Random.seed <- c(57, 8, 4, 7, 43, 3, 0, 9, 24, 13, 53, 3)
  > #
  > set.seed(sum(c(57, 8, 4, 7, 43, 3, 0, 9, 24, 13, 53, 3)))
  > #
  > schpermwage1<-sample(original$chpermwage,size=100000,
  +                      prob=original$wgt/sum(original$wgt),replace=TRUE)
  > schpermwage2<-sample(recent$chpermwage,size=100000,
  +                      prob=recent$wgt/sum(recent$wgt),replace=TRUE)
  > #
  > # wage change (34 wage - 16 wage)
  > #
  > # Try to overlapping densities
  > #
  > nbar <- log(length(schpermwage1), base = 2) + 1
  > kwidth <- diff(range(schpermwage1))/nbar * 0.5
  > swidth <- min(sqrt(var(schpermwage1)),diff(quantile(schpermwage1,c(0.25,0.75)))/1.34)
  > swidth <- 1.06*swidth*length(schpermwage1)^(-1/5)
  > #
  > # Following is the width I used
  > # kwidth <- 1.3*kwidth
  > #
  > # Use Silverman's width
  > #
  > kwidth <- swidth
  > kwidth <- 0.2597
  > dens1 <- density(schpermwage1, n = 500, width=kwidth)
  > plot(x = (dens1$x), y = dens1$y, type = "l", 
  +   xlab = "change in log permanent wage", ylab = "density", 
  +   axes = FALSE,
  +   xlim = c(-1, 3),
  +   ylim=c(0,1.2))
  > title(main="(a)",cex=0.6)
  > axis(side = 1)
  > axis(side = 2)
  > fig1legend <- list(x=c(1.2,1.2),y=c(1.2,1.2))                  
  > legend(fig1legend,lty=1:2,cex=0.5, bty="n",
  + legend=c("original cohort","recent cohort"))       
  > #
  > swidth <- min(sqrt(var(schpermwage2)),diff(quantile(schpermwage2,c(0.25,0.75)))/1.34)
  > swidth <- 1.06*swidth*length(schpermwage2)^(-1/5)
  > #
  > # Following is the width I used
  > # kwidth <- 1.3*kwidth
  > #
  > # Use Silverman's width
  > #
  > kwidth <- swidth
  > kwidth <- 0.455
  > #
  > dens2 <- density(schpermwage2, n = 500, width=kwidth)
  > lines(x = (dens2$x), y = dens2$y, type = "l",lty=2) 
  > #
  > # Now do the Lorenz
  > #
  > swage1 <- sort(recent$chpermwage)
  > swage2 <- sort(original$chpermwage)
  > xout <- (0:1000)/1000
  > #
  > alpha <- seq(along=swage2)/length(swage2)
  > galpha <- cumsum(swage2)/sum(swage2)
  > fn1 <- approx(x=alpha,y=galpha,xout=xout)
  > plot(x = alpha, y = galpha, type = "l", 
  +   xlab = "proportion of population", 
  +   ylab = "proportion of wages", 
  +   ylim=c(0,1.0))
  > legend(x=c(0,0),y=c(1.03,1.03),lty=1:2,cex=0.5, bty="n",
  +  legend=c("original cohort","recent cohort"))
  > #
  > alpha <- seq(along=swage1)/length(swage1)
  > galpha <- cumsum(swage1)/sum(swage1)
  > lines(x = alpha, y = galpha, lty = 2)
  > lines(alpha,alpha,lty=3)
  > title(main="(b)",cex=0.6)
  > #
  > # Now do the RD with smooth=0.4
  > #
  > pdf("Fig2.pdf", width=4.0,height=7.5)
  > 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,
  +         method="bgk",
  +         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,
  +         method="bgk",
  +         ylab="Relative Density",
  +         xlab="Proportion of the Original Cohort")
  > title(main="(b)",cex=0.6)
  > #
  > # Now do the RD with smooth=1.2
  > #
  > pdf("Fig3.pdf", width=4.0,height=7.5)
  > par(err=-1)
  > par(mfrow=c(2,1))
  > #
  > fig3a <- reldist(y=recent$chpermwage,yo=original$chpermwage,ci=FALSE,
  +         smooth=1.2, method="loglik",
  + 	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=1.2,
  + 	yowgt=original$wgt,ywgt=recent$wgt,bar=TRUE,
  + 	yolabs=seq(-1,3,by=0.5),
  + 	ylim=c(0,2.5),cex=0.8,
  +         method="bgk",
  +         ylab="Relative Density",
  +         xlab="Proportion of the Original Cohort")
  > title(main="(b)",cex=0.6)
  > #
  > pdf(file = "Fig4.pdf",width=6.5,height=3.0,paper="USr")
  > 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, smooth=0.4, ci=FALSE,
  + 	ywgt=recent$wgt, yowgt=original$wgt,
  +         method="gam",
  + 	yolabs=seq(-1,3,by=0.5),
  + 	ylim=c(0.5,3.0),
  + 	bar=TRUE, quiet=FALSE,
  + 	xlab="proportion of the original cohort")
  Smoothing using 0.4 
  There were 50 or more warnings (use warnings() to see the first 50)
  > 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="",
  + 	smooth=0.4, ci=FALSE,
  +         method="bgk",
  + 	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)
  Warning message:
  In egvo * log(egv) * dgvo :
    longer object length is not a multiple of shorter object length
  > 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, smooth=0.4, ci=FALSE,
  + 	ywgt=recent$wgt, yowgt=original$wgt,
  + 	show="residual",
  + 	bar=TRUE, quiet=FALSE,
  + 	ylim=c(0.5,3.0), ylab="",
  +         method="bgk",
  + 	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)
  [1] "0.148" "0.183" "0.219" "0.000"
  > format(rpluy(y=recent$chpermwage,yo=original$chpermwage,
  +        ywgt=recent$wgt,yowgt=original$wgt,pvalue=TRUE),
  +        digits=3)
  [1] "1.21e-01" "1.91e-01" "2.62e-01" "5.19e-08"
  > format(rpluy(y=recent$chpermwage,yo=original$chpermwage,
  +        ywgt=recent$wgt,yowgt=original$wgt,upper=TRUE,pvalue=TRUE),
  +        digits=3)
  [1] "1.05e-01" "1.75e-01" "2.46e-01" "6.01e-07"
  > #
  > # 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("Fig5.pdf", width=4.5,height=4.5)
  > 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)),
  +               method="bgk",
  +               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 = "Fig6.pdf",width=6.5,height=3.0,paper="USr")
  > 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, smooth=0.4, ci=FALSE,
  + 	ywgt=recent$wgt, yowgt=original$wgt,
  + 	yolabs=seq(-1,3,by=0.5),
  + 	ylim=c(0.5,3.0),
  +         method="bgk",
  + 	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="",
  + 	smooth=0.4, ci=FALSE,
  +         method="bgk",
  + 	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, smooth=0.4, 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),
  +         method="bgk",
  +         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)
  [1] "0.148" "0.183" "0.219" "0.000"
  > format(rpluy(y=schpermwage1,yo=original$chpermwage,
  +              ywgt=wschpermwage1,yowgt=original$wgt,
  +              pvalue=TRUE),
  +        digits=3)
  [1] "-0.06358" "-0.00729" " 0.04900" " 0.39976"
  > format(rpluy(y=recent$chpermwage,yo=schpermwage1,
  +              ywgt=recent$wgt,yowgt=wschpermwage1,
  +              upper=TRUE, pvalue=TRUE),
  +        digits=3)
  [1] "1.52e-01" "2.04e-01" "2.55e-01" "6.00e-15"
  > #
  > # 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("Fig7.pdf",width=5.5,height=5.5,paper="USr")
  > 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, smooth=0.4,
  + 	ywgt=wgthsr, yowgt=wgthso,
  + 	bar=TRUE, quiet=FALSE,
  + 	ylim=c(0,4),
  +         method="bgk",
  +         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, smooth=0.4,
  + 	ywgt=wgtscr, yowgt=wgtsco,
  + 	bar=TRUE, quiet=FALSE,
  + 	ylim=c(0,4),
  +         method="bgk",
  +         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)
        [,1]     [,2]     [,3]    
   [1,] " 32.56" " 19.85" "-12.71"
   [2,] " 13.59" "  8.87" " -4.72"
   [3,] "  9.82" "  8.77" " -1.05"
   [4,] "  9.53" "  4.97" " -4.56"
   [5,] "  5.46" "  6.91" "  1.45"
   [6,] "  7.32" " 10.02" "  2.71"
   [7,] "  5.96" "  7.50" "  1.55"
   [8,] "  6.22" "  9.09" "  2.87"
   [9,] "  4.35" "  9.51" "  5.15"
  [10,] "  5.19" " 14.50" "  9.31"
  > #
  > collectresults <- c(
  +   wtd.quantile(pwhsr,weight=wgthsr),
  +   wtd.quantile(pwhso,weight=wgthso),
  +   wtd.quantile(pwhsr,weight=wgthsr)-wtd.quantile(pwhso,weight=wgthso),
  +   exp(wtd.quantile(pwhsr,weight=wgthsr)-wtd.quantile(pwhso,weight=wgthso)),
  +   g10hs$entropy,
  +   g10hs$entropy-reldist(y=pwhsr, yo=pwhso, ywgt=wgthsr, yowgt=wgthso, quiet=FALSE,
  +                      method="bgk",
  +                      show="residual",graph=FALSE)$entropy,
  +   reldist(y=pwhsr, yo=pwhso, ywgt=wgthsr, yowgt=wgthso, quiet=FALSE,
  +                      method="bgk",
  +                      show="residual",graph=FALSE)$entropy,
  +   g10hs$rp[2],
  +   g10hs$rpl[2],
  +   g10hs$rpu[2]
  + )
  > #
  > collectresults <- cbind(collectresults,c(
  +   wtd.quantile(pwscr,weight=wgtscr),
  +   wtd.quantile(pwsco,weight=wgtsco),
  +   wtd.quantile(pwscr,weight=wgtscr)-wtd.quantile(pwsco,weight=wgtsco),
  +   exp(wtd.quantile(pwscr,weight=wgtscr)-wtd.quantile(pwsco,weight=wgtsco)),
  +   g10sc$entropy,
  +   g10sc$entropy-reldist(y=pwscr, yo=pwsco, ywgt=wgtscr, yowgt=wgtsco, quiet=FALSE,
  +                      method="bgk",
  +                      show="residual",graph=FALSE)$entropy,
  +   reldist(y=pwscr, yo=pwsco, ywgt=wgtscr, yowgt=wgtsco, quiet=FALSE,
  +                      method="bgk",
  +                      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)
                                 high-school or less some college
  median log recent                        0.6732646  1.138316185
  median log original                      0.9335917  1.128494947
  diff median log                         -0.2603271  0.009821238
  median ratio (recent/original)           0.7707994  1.009869625
  entropy                                  1.9699729  1.216677355
  median entropy                           0.3113690 -0.107023240
  shape entropy                            1.6586038  1.323700595
  MRP                                      0.1538094  0.171613468
  LRP                                      0.2071227  0.270869548
  URP                                      0.1003899  0.073392989
  > binn <- 10
  > pdf(file = "Fig8.pdf",width=6.5,height=3.5,paper="USr")
  > 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.quantile(pwhsr, weight=wgthsr) + wtd.quantile(pwhso, weight=wgthso),
  + yo=pwscr - wtd.quantile(pwscr, weight=wgtscr) + wtd.quantile(pwsco, weight=wgtsco),
  + ywgt=wgthsr, yowgt=wgtscr, binn=binn)
  > #
  > mhsrdhsrscr <- rdeciles(
  +  y=pwhso - wtd.quantile(pwhso, weight=wgthso) + wtd.quantile(pwhsr, weight=wgthsr),
  + yo=pwsco - wtd.quantile(pwsco, weight=wgtsco) + wtd.quantile(pwscr, weight=wgtscr),
  + ywgt=wgthso, yowgt=wgtsco, binn=binn)
  > #
  > m1rdhsrscr <- rdeciles(yo=pwsco, 
  + y=pwhsr - wtd.quantile(pwhsr, weight=wgthsr) + wtd.quantile(pwhso, weight=wgthso),
  + yowgt=wgtsco, ywgt=wgthsr, binn=binn)
  > #
  > m2rdhsrscr <- rdeciles(y=pwhso, 
  + yo=pwscr - wtd.quantile(pwscr, weight=wgtscr) + wtd.quantile(pwsco, weight=wgtsco),
  + ywgt=wgthso, yowgt=wgtscr, binn=binn)
  > #
  > m3rdhsrscr <- rdeciles(y=pwhsr, 
  + yo=pwsco - wtd.quantile(pwsco, weight=wgtsco) + wtd.quantile(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"))
  > 
  > proc.time()
     user  system elapsed 
     9.28    0.68    9.95