library(plyr)
library(readstata13)
library(statar)
library(data.table)
library(plm)
library(ggplot2)
library(lattice)
library(zoo)
library(weights)

# This file computes a set of graphs
# 3(d) for the main paper
# And then Figures D.1. to D.3 and D.5 for the Appendix


rm(list=ls())

min.age <- 25

  binning <- cumsum(c(rep(0.02,5),rep(0.1,8),rep(0.02,5)))
  eps.probs <- cumsum(c(0.025,0.025,0.05,0.4,0.4,0.05,0.025,0.025))
  eps.probs3 <- c(0.025,0.025,0.05,0.4,0.4,0.05,0.025,0.025)

  df <- as.data.table(read.table("Outputs/eta_sample.txt",sep=","))
  eps <- as.data.table(read.table("Outputs/eps_sample.txt",sep=","))
  
  
  names(df) <- c("id","Age","ry","ny")
  names(eps) <- c("id","Age","ry","ny")
  
  
  total.data <- merge(df,eps,by=c("id","Age"))
  total.data$ry <- total.data$ry.x+total.data$ry.y
  total.data$ny <- total.data$ny.x+total.data$ny.y
  
  total.data$ry.x <- NULL
  total.data$ry.y <- NULL
  total.data$ny.x <- NULL
  total.data$ny.y <- NULL
  
  # Variance of log earnings in NL process is variance of total.data
  
  var.nl <- total.data[,var(ry),by="Age"]
  var.nl <- rename(var.nl,c("V1"="varnl"))
  
  # Variance in the data was stored before
  
  var.data <- read.csv("Outputs/earnvarprofile_age_kaplan.csv")
  var.data$X <- NULL
  var.data <- rename(var.data,c("V1"="vardata"))
  
  # Variance in the canonical process is stored by sty_boot.m
  
  var.canon <- read.table("Outputs/variance_earn_canonical.csv",sep=",")
  var.canon$Age <- seq(25,60)
  var.canon <- rename(var.canon,c("V1"="varcanon"))
  
  # Together
  
  all.vars <- merge(var.nl,merge(var.data,var.canon,by="Age"),by="Age")

  all.vars.melt <- melt(all.vars,id="Age")
  
  # Figure 3.d.
  
  source('functions/gon_graphs_export.R')
  folder <- "Outputs/"
  g1 <- grafico.varias.lineas(all.vars.melt,"Age","value","variable","Age","Variance","f_earnvar.pdf")
  
  # Figures D.1. to D.3 and D.5.
  # Note: moments_process_check.m must be run before getting these graphs
  # For the right hand side of figure D.1., the simulation in Matlab must be run
  # before obtaining the graph (commented out)
  
  source('functions/compute_moments.R')