library(lattice)
library(xtable)

load("results.rda")
rm(x)

# Raw material for making tables
tmp <- sort(unique(all$initial.wage))
wages.short <- c(tmp[1], tmp[ceiling(0.5*length(tmp))], tmp[length(tmp)])
tmp <- sort(unique(all$equityfraction))
equityfraction.short <- c(tmp[1], tmp[ceiling(0.5*length(tmp))], tmp[length(tmp)])

# Picture of Pr(poverty)
d <- rbind(
           data.frame(equityfraction=all$equityfraction,
                      initial.wage=all$initial.wage,
                      Pr=all$g0.Pr,
                      g="None"),
           data.frame(equityfraction=all$equityfraction,
                      initial.wage=all$initial.wage,
                      Pr=all$g1.Pr,
                      g="G1 (>Poverty)"),
           data.frame(equityfraction=all$equityfraction,
                      initial.wage=all$initial.wage,
                      Pr=all$g2.Pr,
                      g="G2 (>50% replacement)"),
           data.frame(equityfraction=all$equityfraction,
                      initial.wage=all$initial.wage,
                      Pr=all$g3.Pr,
                      g="G3 (>0% real return)")
           )
table <- NULL
for (w in wages.short) {
    for (e in equityfraction.short) {
        useful <- subset(d, d$equityfraction==e & d$initial.wage==w)$Pr
        table <- rbind(table, c(w, e, useful))
    }
}
colnames(table) <- c("Init. wage", "Eq. frac.", "None",
                     "G1: >Poverty", "G2: >50% repl.",
                     "G3: >0% real")
cat("Probability of poverty table --\n")
table
xtable(table)

pdf("pr_poverty.pdf")
contourplot(Pr ~ equityfraction + initial.wage | g, data=d,
            cuts=7, pretty=TRUE, region = TRUE)

# Picture of median replacement rate
d <- rbind(
           data.frame(equityfraction=all$equityfraction,
                      initial.wage=all$initial.wage,
                      rr.median=all$g0.rr.median,
                      g="None"),
           data.frame(equityfraction=all$equityfraction,
                      initial.wage=all$initial.wage,
                      rr.median=all$g1.rr.median,
                      g="G1 (>Poverty)"),
           data.frame(equityfraction=all$equityfraction,
                      initial.wage=all$initial.wage,
                      rr.median=all$g2.rr.median,
                      g="G2 (>50% replacement)"),
           data.frame(equityfraction=all$equityfraction,
                      initial.wage=all$initial.wage,
                      rr.median=all$g3.rr.median,
                      g="G3 (>0% real return)")
           )
table <- NULL
for (w in wages.short) {
    for (e in equityfraction.short) {
        useful <- subset(d, d$equityfraction==e & d$initial.wage==w)$rr.median
        table <- rbind(table, c(w, e, useful))
    }
}
colnames(table) <- c("Init. wage", "Eq. frac.", "None",
                     "G1: >Poverty", "G2: >50% repl.",
                     "G3: >0% real")
cat("Median replacement rate table --\n")
table
xtable(table)

pdf("median_rr.pdf")
contourplot(rr.median ~ equityfraction + initial.wage | g, data=d,
            cuts=7, pretty=TRUE, region = TRUE)

# Picture of IQR
d <- rbind(
           data.frame(equityfraction=all$equityfraction,
                      initial.wage=all$initial.wage,
                      rr.IQR=all$g0.rr.IQR,
                      g="None"),
           data.frame(equityfraction=all$equityfraction,
                      initial.wage=all$initial.wage,
                      rr.IQR=all$g1.rr.IQR,
                      g="G1 (>Poverty)"),
           data.frame(equityfraction=all$equityfraction,
                      initial.wage=all$initial.wage,
                      rr.IQR=all$g2.rr.IQR,
                      g="G2 (>50% replacement)"),
           data.frame(equityfraction=all$equityfraction,
                      initial.wage=all$initial.wage,
                      rr.IQR=all$g3.rr.IQR,
                      g="G3 (>0% real return)")
           )
table <- NULL
for (w in wages.short) {
    for (e in equityfraction.short) {
        useful <- subset(d, d$equityfraction==e & d$initial.wage==w)$rr.IQR
        table <- rbind(table, c(w, e, useful))
    }
}
colnames(table) <- c("Init. wage", "Eq. frac.", "None",
                     "G1: >Poverty", "G2: >50% repl.",
                     "G3: >0% real")
cat("IQR of replacement rate table --\n")
table
xtable(table)

pdf("rr_IQR.pdf")
contourplot(rr.IQR ~ equityfraction + initial.wage | g, data=d,
            cuts=7, pretty=TRUE, region = TRUE)

# Picture of price
d <- rbind(
           data.frame(equityfraction=all$equityfraction,
                      initial.wage=all$initial.wage,
                      price=all$g0.price,
                      g="None"),
           data.frame(equityfraction=all$equityfraction,
                      initial.wage=all$initial.wage,
                      price=all$g1.price,
                      g="G1 (>Poverty)"),
           data.frame(equityfraction=all$equityfraction,
                      initial.wage=all$initial.wage,
                      price=all$g2.price,
                      g="G2 (>50% replacement)"),
           data.frame(equityfraction=all$equityfraction,
                      initial.wage=all$initial.wage,
                      price=all$g3.price,
                      g="G3 (>0% real return)")
           )
d$price[d$price==0] <- 1
d$logprice <- log(d$price)/log(10)

table <- NULL
for (w in wages.short) {
    for (e in equityfraction.short) {
        useful <- subset(d, d$equityfraction==e & d$initial.wage==w)$logprice
        table <- rbind(table, c(w, e, useful))
    }
}
colnames(table) <- c("Init. wage", "Eq. frac.", "None",
                     "G1: >Poverty", "G2: >50% repl.",
                     "G3: >0% real")
cat("Log price of guarantee table --\n")
table
xtable(table)

pdf("price.pdf")
contourplot(logprice ~ equityfraction + initial.wage | g, data=d,
            cuts=7, pretty=TRUE, region = TRUE)