# bfvariance-figures.ssc
#
#
# Figures for bfvariance paper
#
# Author: Hans-Peter Kohler
# email:  kohler@demogr.mpg.de


# calculate standardized fertility rates

# "gd" and "nfx" are just interim objects used by the plot functions
gd.1 <- se.os.res[[1]]
gd.2 <- se.os.res[[2]]
gd.3 <- se.os.res[[3]]
gd.a <- se.all.res


nfx <- se.nfx.1
nfx.std <- se.nfx.1
p.age <- se.age + 0.5


### Figure 1 in KP
p.cex <- 0.4
leg.cex <- 0.53
c.lwd <- 1.5
par(mfrow=c(2,3))
#observed data
t.r.bf <- c(NA, 0.5 * (gd.1$mab[3:length(gd.1$years)] - 
                       gd.1$mab[1:(length(gd.1$years)-2)]), NA)
matplot(1900+gd.1$years, cbind(gd.1$tfr,gd.1$tfr.s,gd.1$tfr.adj.bf), type="n",
        lwd=c.lwd,
	col=1, lty=c(2,1,3),
	xlab="year",
	ylab="TFR")
legend(1975,1.05,c("observed TFR", "obs. TFR, smoothed", "adjusted TFR (BF)"), 
	lty=c(2,1,3),lwd=c.lwd,
	bty="n", cex=leg.cex)
matlines(1900+gd.1$years, cbind(gd.1$tfr,gd.1$tfr.s,gd.1$tfr.adj.bf), lwd=c.lwd,
	col=1, lty=c(2,1,3))
title(main="observed TFR and adjusted TFR (BF), order 1", cex=p.cex)

bblev <- .66
ttlev <- .675
lines(c(1985,1990), c(bblev,bblev),type="b",pch="|",lty=1)
text(1987.5,ttlev,"'boom'", cex=leg.cex, adj=.5)
lines(c(1990,1995), c(bblev,bblev),type="b",pch="|",lty=1)
text(1992.5,ttlev,"'bust'", cex=leg.cex, adj=.5)
rm(bblev,ttlev)

matplot(1900+gd.1$years, cbind(gd.1$mab, gd.1$mab.s), type="n",lwd=c.lwd,
        col=1,lty=c(2,1),
	xlab="year",
	ylab="mu")
title(main="observed mean age at birth (mu), order 1", cex=p.cex)
legend(1975,27,c("observed MAFB", "obs. MAFB, smoothed"), 
       lty=c(2,1),lwd=c.lwd,
       bty="n", cex=leg.cex)
matlines(1900+gd.1$years, cbind(gd.1$mab, gd.1$mab.s), lwd=c.lwd,
         col=1,lty=c(2,1))

bblev <- 24.5
ttlev <- 24.6
lines(c(1985,1990), c(bblev,bblev),type="b",pch="|",lty=1)
text(1987.5,ttlev,"'boom'", cex=leg.cex, adj=.5)
lines(c(1990,1995), c(bblev,bblev),type="b",pch="|",lty=1)
text(1992.5,ttlev,"'bust'", cex=leg.cex, adj=.5)
rm(bblev,ttlev)


matplot(1900+gd.1$years, cbind(sqrt(gd.1$vab),sqrt(gd.1$vab.s)), type="n",lwd=c.lwd,
	col=1,lty=c(2,1),
	xlab="year",
	ylab="std. deviation")
title(main="obs. std. deviation of age at birth, order 1", cex=p.cex)
legend(1975,4.9,c("observed std. deviation", "std. deviation, smoothed "), 
       lty=c(2,1),lwd=c.lwd,
       bty="n", cex=leg.cex)
matlines(1900+gd.1$years, cbind(sqrt(gd.1$vab),sqrt(gd.1$vab.s)), lwd=c.lwd,
         col=1,lty=c(2,1))
bblev <- 4.48
ttlev <- 4.5
lines(c(1985,1990), c(bblev,bblev),type="b",pch="|",lty=1)
text(1987.5,ttlev,"'boom'", cex=leg.cex, adj=.5)
lines(c(1990,1995), c(bblev,bblev),type="b",pch="|",lty=1)
text(1992.5,ttlev,"'bust'", cex=leg.cex, adj=.5)
rm(bblev,ttlev)



# plot fertility schedules in bottom panel of Figure 1
yr.min <- 84
yr.max <- 96
plot.yr <- gd.1$years >= yr.min & gd.1$years <= yr.max & gd.1$years %% 3 == 0
n.y <- sum(plot.yr)

#observed nfx
matplot(p.age,nfx[,plot.yr],
	type="l", lwd=c.lwd, col=1, lty=1:n.y,
	xlab="age",ylab="nFx")
title(main="observed age-specific fertility rates nFx", cex=p.cex)
legend(35,0.07, paste(1900+gd.1$years[plot.yr]),lty=1:n.y, lwd=c.lwd,
       col=1, cex=leg.cex, bty="n")

# standarize fertit schedule to same quantum
for (i in 1:dim(nfx.std)[1]){
  nfx.std[i,] <-  nfx.std[i,]/gd.1$tfr
}

#mean-standardize nfx
age.mstd <- matrix(p.age,length(p.age),n.y,byrow=F) - 
  matrix(gd.1$mab[plot.yr],length(p.age),n.y,byrow=T) +
  gd.1$mab[gd.1$years == yr.min]

matplot(age.mstd,nfx.std[,plot.yr],
	type="l", col=1, lty=1:n.y, lwd=c.lwd,	
	xlab="age",ylab="nFx")
title(main="mean- and level-standardized nFx", cex=p.cex)
legend(35,0.08, paste(1900+gd.1$years[plot.yr]),lty=1:n.y,lwd=c.lwd,
       col=1, cex=leg.cex, bty="n")

#mean and variance standardize
v.ratio <- sqrt(gd.1$vab[plot.yr]/gd.1$vab[gd.1$years==yr.min])
age.vstd <- ((matrix(p.age,length(p.age),n.y,byrow=F) - 
              matrix(gd.1$mab[plot.yr],length(p.age),n.y,byrow=T)) /
             matrix(v.ratio,length(p.age),n.y,byrow=T))+
              gd.1$mab[gd.1$years == yr.min]
nfx.vstd <- nfx.std[,plot.yr]*matrix(v.ratio,length(p.age),n.y,byrow=T)
matplot(age.vstd,nfx.vstd,
	type="l", col=1, lty=1:n.y,lwd=c.lwd,
	xlab="age",ylab="nFx")
title(main="mean-, level-, and variance-standardized nFx", cex=p.cex)
legend(35,0.08, paste(1900+gd.1$years[plot.yr]),lty=1:n.y,lwd=c.lwd,
	col=1, cex=leg.cex, bty="n")


# Figure 3 in KP
yr.min <- 75
c.lwd <- 1.8
par(mfrow=c(2,3))

#plot observed variance effect
plot(1900+gd.1$years,gd.1$delta, type="l", lwd=c.lwd, col=1,
     xlab="year",xlim=c(1975,1995),
     ylab="change in std. deviation")
lines(c(1900,2000),c(0,0), lty=2)
legend(1975,0.012,c("delta"), lty=1,lwd=c.lwd,bty="n", cex=leg.cex)
title(main="annual change in std. deviation (order 1)", cex=p.cex)

bblev <- -.0012
ttlev <- -.0006
lines(c(1985,1990), c(bblev,bblev),type="b",pch="|",lty=1)
text(1987.5,ttlev,"'boom'", cex=leg.cex, adj=.5)
lines(c(1990,1995), c(bblev,bblev),type="b",pch="|",lty=1)
text(1992.5,ttlev,"'bust'", cex=leg.cex, adj=.5)
rm(bblev,ttlev)


matplot(1900+gd.1$years, cbind(gd.1$r.bf,gd.1$r.bf.biascorr,gd.1$gamma), type="l",
        lwd=c.lwd,
	col=1, lty=c(3,2,1), xlim=1900+c(yr.min,yr.max), ylim=c(0,0.255),
	xlab="year",
	ylab="tempo")
legend(1975,0.255,c("gamma (with variance effect)",
                    "r (Bongaarts-Feeney)", "r (BF), bias corrected"),
       lty=c(1,3,2),lwd=c.lwd,
       bty="n", cex=leg.cex)
title(main="annual change in mean age (order 1)", cex=p.cex)
bblev <- 0
ttlev <- .008
lines(c(1985,1990), c(bblev,bblev),type="b",pch="|",lty=1)
text(1987.5,ttlev,"'boom'", cex=leg.cex, adj=.5)
lines(c(1990,1995), c(bblev,bblev),type="b",pch="|",lty=1)
text(1992.5,ttlev,"'bust'", cex=leg.cex, adj=.5)
rm(bblev,ttlev)


#just an empty plot
plot(c(0,1),c(0,1),type="n",ylab="")

matplot(1900+gd.1$years, cbind(gd.1$mab.s, gd.1$mab.adj.var), type="l",lwd=c.lwd,
        col=1,lty=c(1,2),xlim=c(1975,1995),
	xlab="year",
	ylab="mu")
title(main="mean age at birth", cex=p.cex)
legend(1975,27,c("observed MAFB", "corrected MAFB"), 
       lty=c(1,2),lwd=c.lwd,
       bty="n", cex=leg.cex)

arrows(1987.5,25.1,1987.5,26.1, open=T, size=.1)
text(1987.5,25,"'plateau' in MAB",cex=leg.cex)

bblev <- 24.5
ttlev <- 24.63
lines(c(1985,1990), c(bblev,bblev),type="b",pch="|",lty=1)
text(1987.5,ttlev,"'boom'", cex=leg.cex, adj=.5)
lines(c(1990,1995), c(bblev,bblev),type="b",pch="|",lty=1)
text(1992.5,ttlev,"'bust'", cex=leg.cex, adj=.5)
rm(bblev,ttlev)


matplot(1900+gd.1$years, cbind(gd.1$tfr.adj.var,gd.1$tfr.adj.bf,gd.1$tfr.s),
	type="n",lwd=1, xlim=1900+c(yr.min,yr.max),lwd=c.lwd,
	col=1, lty=c(1,3,2),
	xlab="year",
	ylab="TFR")
legend(1900+yr.min,1.07,c("adj. TFR w. var. effect",
                          "adj. TFR (BF)","observed TFR"), 
       lwd=c.lwd,lty=c(1,3,2), bty="n", cex=leg.cex)
matlines(1900+gd.1$years, cbind(gd.1$tfr.adj.var,gd.1$tfr.adj.bf,gd.1$tfr.s),
         lwd=c.lwd,, col=1, lty=c(1,3,2))
title(main="observed and adjusted TFR (order 1)", cex=p.cex)
bblev <- .66
ttlev <- .677
lines(c(1985,1990), c(bblev,bblev),type="b",pch="|",lty=1)
text(1987.5,ttlev,"'boom'", cex=leg.cex, adj=.5)
lines(c(1990,1995), c(bblev,bblev),type="b",pch="|",lty=1)
text(1992.5,ttlev,"'bust'", cex=leg.cex, adj=.5)
rm(bblev,ttlev)



#just an empty plot
plot(c(0,1),c(0,1),type="n", ylab="")



# Figure 4 in KP

par(mfrow=c(2,3))
matplot(1900+gd.a$years, cbind(gd.1$r.bf,gd.2$r.bf,gd.3$r.bf), 
	type="n",lwd=c.lwd,
	xlim=c(1900+yr.min, 1900+yr.max),ylim=c(-0.03,0.2),
	xlab="year",
	ylab="annual change in MAB")
matlines(1900+gd.a$years, cbind(gd.1$r.bf,gd.2$r.bf,gd.3$r.bf), 
         type="l", lwd=c.lwd,
         col=1, lty=c(1,2,3))
title(main="annual change in MAB (BF-method)", cex=p.cex)
legend(1975,0.02,c("first births", "second births",
                   "third births"), lty=c(1,2,3), lwd=c.lwd,
       bty="n", cex=p.cex)

matplot(1900+gd.a$years, cbind(gd.1$gamma,gd.2$gamma,gd.3$gamma), 
	type="n",lwd=c.lwd,
	xlim=c(1900+yr.min, 1900+yr.max), ylim=c(-0.03,0.2),
	xlab="year",
	ylab="annual change in MAB")
matlines(1900+gd.a$years, cbind(gd.1$gamma,gd.2$gamma,gd.3$gamma), 
         type="l", lwd=c.lwd,
         col=1, lty=c(1,2,3))
title(main="annual change in MAB w. variance effect", cex=p.cex)
legend(1975,0.02,c("first births", "second births",
                   "third births"), lty=c(1,2,3), lwd=c.lwd,
       bty="n", cex=p.cex)

text(1989,0.19,"policy effect 1", cex=p.cex)
text(1981.5,0.035,"policy effect 2", cex=p.cex)
arrows(c(1987.75,1988.25),c(0.18,.18),c(1980.9,1983.2), c(.135,.085), open=T, size=.1)
arrows(1985,0.034,1988,0.034, open=T, size=.1)