# Performance of estimation

source('functions/gon_graphs_export.R')

source('functions/moment_functions.R')

setwd('Outputs/')

folder <- ""

min.age <- 25

# First, perform some computations in the SIMULATED DATA after applying
# the ABB decomposition, and compute the associated moments.

td <- subset(total.data,select=c("id","Age","ry","ny"))
td$ng <- td$ny - td$ry

setkeyv(td,c("Age","id"))
eje <- td[,xtile(ry,probs=binning),by=Age]
td$pc <- eje$V1


td$pc <- mapvalues(td$pc,from=seq(1,length(binning)),
                   to=round(binning*100))

menv <- new.env()
uncond.moments(td,menv,min.age)
td <- gen.agroups(td,2)


td <- subset(td,ry>-6) # trim very large outliers
                      # (simulation-dependent)
td <- subset(td,ry<5)

#Computing condmoments

cond.moments(td,menv)

# Do the same for the PERSISTENT COMPONENT ONLY

setkeyv(df,c("Age","id"))
eje <- df[,xtile(ry,probs=binning),by=Age]
df$pc <- eje$V1
df$pc <- mapvalues(df$pc,from=seq(1,length(binning)),
                   to=round(binning*100))

da <- df

penv <- new.env()
uncond.moments(da,penv,min.age) 
da <- gen.agroups(da,2)
da$ng <- da$ny - da$ry

cond.moments(da,penv)


tenv <- new.env() #transitory component now
uncond.moments(eps,tenv,min.age)


# Load from moments_process_check.m the implied moments for the discretized processes


sd.bis <- read.csv('mm_sd.csv',header=F)
sd.bis$pc <- binning*100
sd.bis <- melt(sd.bis,id="pc")
sd.bis$variable <- mapvalues(sd.bis$variable,from=c(paste("V",seq(1,3),sep="")),to=
                               c("25-34","35-44","45-54"))


sk.bis <- read.csv('mm_sk.csv',header=F)
sk.bis$pc <- binning*100
sk.bis <- melt(sk.bis,id="pc")
sk.bis$variable <- mapvalues(sk.bis$variable,from=c(paste("V",seq(1,3),sep="")),to=
                               c("25-34","35-44","45-54"))


kt.bis <- read.csv('mm_kt.csv',header=F)
kt.bis$pc <- binning*100
kt.bis <- melt(kt.bis,id="pc")
kt.bis$variable <- mapvalues(kt.bis$variable,from=c(paste("V",seq(1,3),sep="")),to=
                               c("25-34","35-44","45-54"))


beps2 <- tenv$usd

sd.tra <-  read.csv('mm_sdtrans.csv',header=F)
beps2$Discretized <- sd.tra$V1
beps2 <- rename(beps2,c("V1"="Data"))

ceps2 <- tenv$uskew

sk.tra <-  read.csv('mm_sktrans.csv',header=F)
ceps2$Discretized <- sk.tra$V1
ceps2 <- rename(ceps2,c("V1"="data"))

deps2 <- tenv$ukurt

kt.tra <-  read.csv('mm_kttrans.csv',header=F)
deps2$Discretized <- kt.tra$V1
deps2 <- rename(deps2,c("V1"="data"))

beps3 <- tenv$usd

sd.tra3 <-  read.csv('mm_sdtrans2.csv',header=F)
beps3$Discretized <- sd.tra3$V1
beps3 <- rename(beps3,c("V1"="Data"))

ceps3 <- tenv$uskew

sk.tra3 <-  read.csv('mm_sktrans2.csv',header=F)
ceps3$Discretized <- sk.tra3$V1
ceps3 <- rename(ceps3,c("V1"="data"))

deps3 <- tenv$ukurt

kt.tra3 <-  read.csv('mm_kttrans2.csv',header=F)
deps3$Discretized <- kt.tra3$V1
deps3 <- rename(deps3,c("V1"="data"))




# Figure D.1., central panel. The left-hand side panel is in Figure 1,
# and the right-hand side panel can be replicated from the simulation in the Model
# section

d1b <- cond.age.moments(menv$csd,"pc","condsd","agroup","Percentile of previous earnings",
                             "Std dev of log earnings change","f_sumsim_sd.pdf",c(0,1),7,NULL,26,NULL,1.5)

d1e <- cond.age.moments(menv$cskew,"pc","V1","agroup","Percentile of previous earnings",
                             "Skewness of log earnings change","f_sumsim_sk.pdf",c(-4,1),7,NULL,26,NULL,1.5)

d1h <- cond.age.moments(menv$ckurt,"pc","V1","agroup","Percentile of previous earnings",
                             "Kurtosis of log earnings change","f_sumsim_kt.pdf",c(0,30),7,NULL,26,NULL,1.5)


# Figure D.2.

d2a <- cond.age.moments(penv$csd,"pc","condsd","agroup","Percentile of previous earnings",
                                 "Standard deviation","condsd.pdf",c(0.17,0.63),7,NULL,26,NULL,1.5)

d2b <- cond.age.moments(penv$cskew,"pc","V1","agroup","Percentile of previous earnings",
                                  "Conditional skewness","condskew.pdf",NULL,7,NULL,26,NULL,1.5)

d2c <- cond.age.moments(penv$ckurt,"pc","V1","agroup","Percentile of previous earnings",
                                  "Conditional kurtosis","condkurt.pdf",c(0,30),7,NULL,26,NULL,1.5)


d2d <- cond.age.moments(sd.bis,"pc","value","variable","Percentile of previous earnings",
                                 "Standard deviation","condsd_discr.pdf",c(0.17,0.63),7,NULL,26,NULL,1.5)

d2e <- cond.age.moments(sk.bis,"pc","value","variable","Percentile of previous earnings",
                                  "Conditional skewness","condskew_discr.pdf",NULL,7,NULL,26,NULL,1.5)

kt.bis <- rename(kt.bis,c("variable"="Age group"))
d2f <- cond.age.moments(kt.bis,"pc","value","Age group","Percentile of previous earnings",
                                  "Conditional kurtosis","condkurt_discr.pdf",c(0,30),8,NULL,27,2.1,1.5)

# Figure D.3., left panel. The right panel can be obtained by running moments_process_check.m
# on the corresponding state spaces + transition matrices

beps2 <- melt(beps2,id="Age")
beps2$variable <- mapvalues(beps2$variable,from=c("Data","Discretized"),to=c("Transitory earnings",
                                                                             "Discretized tr. earnings"))
d3a <- grafico.varias.lineas.bn(beps2,"Age","value","variable","Age",
                               "Standard deviation","f_mm_sdtrans.pdf",c(0,0.5),NULL,0.3,20)

ceps2 <- melt(ceps2,id="Age")
ceps2$variable <- mapvalues(ceps2$variable,from=c("data","Discretized"),to=c("Transitory earnings",
                                                                             "Discretized tr. earnings"))
d3c <- grafico.varias.lineas.bn(ceps2,"Age","value","variable","Age","Skewness","f_mm_sktrans.pdf",
                               c(-4,0),NULL,0.3,20)

deps2 <- melt(deps2,id="Age")
deps2$variable <- mapvalues(deps2$variable,from=c("data","Discretized"),to=c("Transitory earnings",
                                                                             "Discretized tr. earnings"))
d3e <- grafico.varias.lineas.bn(deps2,"Age","value","variable","Age","Kurtosis","f_mm_kttrans.pdf",
                               c(0,65),NULL,0.3,20)


beps3 <- melt(beps3,id="Age")
beps3$variable <- mapvalues(beps3$variable,from=c("Data","Discretized"),to=c("Transitory earnings",
                                                                             "Discretized tr. earnings"))
d3b <- grafico.varias.lineas.bn(beps3,"Age","value","variable","Age",
                                "Standard deviation","f_mm_sdtrans2.pdf",c(0,0.5),NULL,0.3,20)

ceps3 <- melt(ceps3,id="Age")
ceps3$variable <- mapvalues(ceps3$variable,from=c("data","Discretized"),to=c("Transitory earnings",
                                                                             "Discretized tr. earnings"))
d3d <- grafico.varias.lineas.bn(ceps3,"Age","value","variable","Age","Skewness","f_mm_sktrans2.pdf",
                                c(-4,0),NULL,0.3,20)

deps3 <- melt(deps3,id="Age")
deps3$variable <- mapvalues(deps3$variable,from=c("data","Discretized"),to=c("Transitory earnings",
                                                                             "Discretized tr. earnings"))
d3f <- grafico.varias.lineas.bn(deps3,"Age","value","variable","Age","Kurtosis","f_mm_kttrans2.pdf",
                                c(0,65),NULL,0.3,20)

# Figure D.5., left panels. The right hand side panels can be obtained by running moments_process_check.m
# on the corresponding state spaces + transition matrices

sd.unc <-  read.csv('mm_sdunc.csv',header=F)
bb2 <- penv$usd
bb2$Discretized <- sd.unc$V1
bb2 <- rename(bb2,c("V1"="Data"))

sk.unc <-  read.csv('mm_skunc.csv',header=F)
cc2 <- penv$uskew
cc2$Discretized <- sk.unc$V1
cc2 <- rename(cc2,c("V1"="Data"))

kt.unc <-  read.csv('mm_ktunc.csv',header=F)
dd2 <- penv$ukurt
dd2$Discretized <- kt.unc$V1
dd2 <- rename(dd2,c("V1"="Data"))

sd.unc3 <-  read.csv('mm_sdunc2.csv',header=F)
bb3 <- penv$usd
bb3$Discretized <- sd.unc3$V1
bb3 <- rename(bb3,c("V1"="Data"))

sk.unc3 <-  read.csv('mm_skunc2.csv',header=F)
cc3 <- penv$uskew
cc3$Discretized <- sk.unc3$V1
cc3 <- rename(cc3,c("V1"="Data"))

kt.unc3 <-  read.csv('mm_ktunc2.csv',header=F)
dd3 <- penv$ukurt
dd3$Discretized <- kt.unc3$V1
dd3 <- rename(dd3,c("V1"="Data"))

bbgraph <- melt(bb2,id="Age")
bbgraph$variable <- mapvalues(bbgraph$variable,from=c("Data","Discretized"),to=c("Persistent earnings"
                                                                                 ,"Discretized pers. earnings"))
d5a <- grafico.varias.lineas.bn(bbgraph,"Age","value","variable","Age","Standard deviation","f_mm_sdunc.pdf",
                               c(0,1))

ccgraph <- melt(cc2,id="Age")
ccgraph$variable <- mapvalues(ccgraph$variable,from=c("Data","Discretized"),to=c("Persistent earnings"
                                                                                 ,"Discretized pers. earnings"))
d5c <- grafico.varias.lineas.bn(ccgraph,"Age","value","variable","Age","Skewness","f_mm_skunc.pdf",
                               c(-1.25,0))

ddgraph <- melt(dd2,id="Age")
ddgraph$variable <- mapvalues(ddgraph$variable,from=c("Data","Discretized"),to=c("Persistent earnings"
                                                                                 ,"Discretized pers.earnings"))
d5e <- grafico.varias.lineas.bn(ddgraph,"Age","value","variable","Age","Kurtosis","f_mm_ktunc.pdf",
                               c(3,7))

bbgraph2 <- melt(bb3,id="Age")
bbgraph2$variable <- mapvalues(bbgraph2$variable,from=c("Data","Discretized"),to=c("Persistent earnings"
                                                                                 ,"Discretized pers. earnings"))
d5b <- grafico.varias.lineas.bn(bbgraph2,"Age","value","variable","Age","Standard deviation","f_mm_sdunc2.pdf",
                                c(0,1))

ccgraph2 <- melt(cc3,id="Age")
ccgraph2$variable <- mapvalues(ccgraph2$variable,from=c("Data","Discretized"),to=c("Persistent earnings"
                                                                                 ,"Discretized pers. earnings"))
d5d <- grafico.varias.lineas.bn(ccgraph2,"Age","value","variable","Age","Skewness","f_mm_skunc2.pdf",
                                c(-1.25,0))

ddgraph2 <- melt(dd3,id="Age")
ddgraph2$variable <- mapvalues(ddgraph2$variable,from=c("Data","Discretized"),to=c("Persistent earnings"
                                                                                 ,"Discretized pers.earnings"))
d5f <- grafico.varias.lineas.bn(ddgraph2,"Age","value","variable","Age","Kurtosis","f_mm_ktunc2.pdf",
                                c(3,7))


# Figure D.1., discretized components, requires specific treatment because persistent+
# transitory component need to be added
# Take it from simulation in Matlab
# Commented out here to get the code running even before simulation was run

# 
# ## Now moments for ALL earnings
# 
# 
# file.name <- paste("../../Model/Outputs/dfp_simulation_choose128_bc2.csv",sep="")
# 
# dat <- as.data.table(read.table(file.name,sep=","))
# dat <- rename(dat,c("V1"="id","V2"="Age","V3"="ry","V4"="rye","V5"="sav","V6"="cons","V7"="newlogy"))
# 
# dat$ry <- dat$newlogy
# dat[, nexty:=shift(ry,type="lead"),by="id"]
# dat$toty <- dat$ry + log(dat$rye)
# dat[, nextoty:=shift(toty,type="lead"),by="id"]
# dat$ng <- dat$nextoty - dat$toty
# setkeyv(dat,c("Age","id"))
# eje <- dat[,xtile(toty,probs=binning),by=Age]
# dat$pc <- eje$V1
# dat$pc <- mapvalues(dat$pc,from=seq(1,length(binning)),
#                    to=round(binning*100))
# dat <- subset(dat,is.na(nextoty)==0)
# 
# dat$ny <- NULL
# dat$ry <- NULL
# dat <- rename(dat,c("toty"="ry","nextoty"="ny"))
# senv <- new.env()
# dat <- gen.agroups(dat,2)
# uncond.moments(dat,senv,min.age)
# cond.moments(dat,senv)
# 
# discr.sd <- cond.age.moments(senv$csd,"pc","condsd","agroup","Percentile of previous earnings",
#                              "Std dev of log earnings change","f_discr_sd.pdf",c(0,1),8,NULL,26,64,1.5)
# 
# discr.sk <- cond.age.moments(senv$cskew,"pc","V1","agroup","Percentile of previous earnings",
#                              "Skewness of log earnings change","f_discr_sk.pdf",c(-4,1),8,NULL,26,12.5,1.5)
# 
# discr.kt <- cond.age.moments(senv$ckurt,"pc","V1","agroup","Percentile of previous earnings",
#                              "Kurtosis of log earnings change","f_discr_kt.pdf",c(0,30),8,NULL,26,2.1,1.5)
