Risk_pensionRisk(5).R

# This script is for modeling tax revenue of stylized governments

# (4) workforce +0.5%; O10d replaces O30p

#**********************************************************************
#                           Packages                               ####
#**********************************************************************


# library(plyr)

library(broom)
library(readxl)
library(magrittr)
library(ggrepel)
library(stringr)
library(forcats)
library(grid)
library(gridExtra)
library(scales)
library(knitr)

library(xlsx)

# packages for econometric and time series modeling
# library(plm)
library(astsa)    # companion package
library(TSA)      # companion package;  arimax: flexible transfer function model
library(tseries)  #
library(forecast) # Arima
library(MSwM)
library(TTR)
library(dynlm)
library(broom)

#library(MSBVAR)

# packages for ts
library(zoo)
library(xts)


library(timetk)
library(tidyquant)

library(lubridate)
library(feather)

library(psych) # describe
library(tidyverse)

options(tibble.print_max = 60, tibble.print_min = 60)


# check tidyquant, timetk, sweep (broom ), tibbletime
# Intro to zoo  cran.r-project.org/web/packages/zoo/vignettes/zoo-quickref.pdf
# sweep: http://www.business-science.io/code-tools/2017/07/09/sweep-0-1-0.html

#**********************************************************************
#                     Global settings and tools                    ####
#**********************************************************************
dir_data_raw <- "data_raw/"
dir_data_out <- "data_out/"
dir_fig_out <- "outputs_report/"
IO_folder <- "C:/Git/PenSim-Projects/Model_Main/IO_penSimMacro"

# NBER recession periods, post-WWII
recessionPeriods <- 
	matrix(c(
		1953+2/4, 1954+2/4,
		1957+3/4, 1958+2/4,
		1960+2/4, 1961+1/4,
		1969+4/4, 1970+4/4,
		1973+4/4, 1975+1/4,
		1980+1/4, 1980+3/4,
		1981+3/4, 1982+4/4,
		1990+3/4, 1991+1/4,
		2001+1/4, 2001+4/4,
		2007+4/4, 2009+2/4
	) , ncol = 2, byrow = T) %>% 
	as.data.frame() %>% 
	rename(peak =   V1,
				 trough = V2) %>% 
	mutate(peak = peak - 1/4,
				 trough = trough - 1/4)


get_logReturn <- function(x){
	if(any(x <= 0, na.rm = TRUE)) stop("Nagative value(s)")
	log(x/dplyr::lag(x))
}




# RIG colors and theme
RIG.blue  <- "#003598"
RIG.red   <- "#A50021"
RIG.green <- "#009900"
RIG.yellow <- "#FFFF66"
RIG.purple <- "#9966FF"
RIG.yellow.dark <- "#ffc829"
RIG.orange <- "#fc9272"

demo.color6 <- c(RIG.red,
								 RIG.orange,
								 RIG.purple,
								 RIG.green ,
								 RIG.blue,
								 RIG.yellow.dark)


RIG.theme <- function() {
	theme(
		panel.grid.major.x = element_blank(),
		panel.grid.minor.x = element_blank(),
		panel.grid.minor.y = element_blank(),
		panel.grid.major.y = element_line(size = 0.5, color = "gray80"),
		plot.title = element_text(hjust = 0.5),
		plot.subtitle = element_text(hjust = 0.5),
		plot.caption = element_text(hjust = 0, size = 9)
	)
}

RIG.themeLite <- function() {
	theme(
		plot.title = element_text(hjust = 0.5),
		plot.subtitle = element_text(hjust = 0.5),
		plot.caption = element_text(hjust = 0, size = 9)
	)
}


get_geoReturn <- function(x) prod(1 + x)^(1/length(x)) - 1


color_PIT <- RIG.green
color_salesgen <- "blue"
color_salessel <- "deepskyblue2"
color_other <- RIG.purple
color_propertyLoc <- RIG.yellow.dark
color_GDP <- "gray50"



#**********************************************************************
#                          Outline                                 ####
#**********************************************************************



#**********************************************************************
#                         Loading data                             ####
#**********************************************************************

# Economic and revenue data: simulation with forward-looking assumptions
load(paste0(dir_data_out, "GovFin_sim_forward_tot.RData"))

# Outputs of pension Finance  
IO_folder

get_results <- function(IO_folder, Pattern = "^Outputs"){
	
	fn <- function(x) {
		load(paste0(IO_folder, "/", x))
		
		# if("results.t7" %in% names(outputs_list)){
		#   df_out <- bind_rows(outputs_list$results,
		#                       outputs_list$results.t7,
		#                       outputs_list$results.xt7)
		#   return(df_out)
		# } else {
		#   return(outputs_list$results)
		# }
		
		return(outputs_list$results)
		
	}
	
	file_select <- dir(IO_folder, Pattern)
	results_all <- plyr::adply(file_select, 1, fn) %>% select(-X1)
}

results_all <- get_results(IO_folder) %>% select(runname, sim, year, everything())
#results_all$runname %>% unique
results_all
#**********************************************************************
#                         Organizing data                          ####
#**********************************************************************

# Assumption on ERC as a percent of government tax revenue 
ERC_taxRev_init <- 0.05

infl <- 0.02

df_sim %>% head




df_taxRev_real <- 
df_sim %>% 
	select(sim, year, 
				 taxLevelReal_PITState_tot_a1,
				 #taxLevelReal_PITState_tot_a2,
				 
				 taxLevelReal_salesState_tot_a1,
				 #taxLevelReal_salesState_tot_a2,
				 
				 taxLevelReal_local_tot_a1,
				 #taxLevelReal_local_tot_a2,
				 
				 taxLevelReal_GDPState_tot_a1
				 
				 ) %>% 
	mutate_at(
		vars(
			taxLevelReal_PITState_tot_a1,
			#taxLevelReal_PITState_tot_a2,
			
			taxLevelReal_salesState_tot_a1,
			#taxLevelReal_salesState_tot_a2,
			
			taxLevelReal_local_tot_a1,
			#taxLevelReal_local_tot_a2,
			
			taxLevelReal_GDPState_tot_a1
			),
		
		funs(ifelse(year == 1, 1, dplyr::lag(.)))
	)
	
	

df_pension <- 
	results_all %>% 
	  filter(sim>0) %>% 
		select(runname, sim, year, C, ERC)

# assuming the ERC under slow repayment of UAAL and 7.5% discount rate is 3% of total tax revenue in year 1
v <- df_pension %>% filter(year == 1, sim %in% 1)
v
taxRev_init <- v[v$runname == "A_O30p_port70_30","ERC"]/ERC_taxRev_init

v %>% mutate(pct = 100* ERC / taxRev_init)





df <- 
	df_pension %>% 
	left_join(df_taxRev_real) %>% 
	group_by(runname, sim) %>% 
	mutate(# Contributions in real terms
		     price_index = (1+infl)^(year-1),
				 C_real   = C / price_index,
				 ERC_real = ERC/ price_index
				 ) %>% 
	mutate_at(
			vars(# Tax revenue 
				taxLevelReal_PITState_tot_a1,
				#taxLevelReal_PITState_tot_a2,
				
				taxLevelReal_salesState_tot_a1,
				#taxLevelReal_salesState_tot_a2,
				
				taxLevelReal_local_tot_a1,
				#taxLevelReal_local_tot_a2,
				
				taxLevelReal_GDPState_tot_a1
				
				),
			
			funs(. * taxRev_init)
	) %>% 
	mutate(
		ERC_tax_PITState_a1   = ERC_real / taxLevelReal_PITState_tot_a1,
		#ERC_tax_PITState_a2   = ERC_real / taxLevelReal_PITState_tot_a2,
		
		ERC_tax_salesState_a1 = ERC_real / taxLevelReal_salesState_tot_a1,
		#ERC_tax_salesState_a1 = ERC_real / taxLevelReal_salesState_tot_a2,
		
		ERC_tax_local_a1      = ERC_real / taxLevelReal_local_tot_a1,
		#ERC_tax_local_a1      = ERC_real / taxLevelReal_local_tot_a2,
		
		ERC_tax_GDPState_a1   = ERC_real / taxLevelReal_GDPState_tot_a1
	)

# x <- df %>% filter(sim %in% 1:3, year == 1)

na2zero <- function(x){replace(x, is.na(x), 0)}


df_risk <- 
	df %>%
	group_by(runname, sim) %>%
	mutate(
	
		ERC_tax_PITState_a1_high1   = cumany(ERC_tax_PITState_a1 >= ERC_tax_PITState_a1[year == 1] + 0.05 ),
		ERC_tax_PITState_a1_high2   = cumany(ERC_tax_PITState_a1 >= 0.12),
		ERC_tax_PITState_a1_hike    = cumany(na2zero(ERC_tax_PITState_a1 - lag(ERC_tax_PITState_a1, 2) >= 0.03)),

		ERC_tax_salesState_a1_high1 = cumany(ERC_tax_salesState_a1>= ERC_tax_salesState_a1[year == 1] + 0.05),
		ERC_tax_salesState_a1_high2 = cumany(ERC_tax_salesState_a1>= 0.12),
		ERC_tax_salesState_a1_hike  = cumany(na2zero(ERC_tax_salesState_a1 - lag(ERC_tax_salesState_a1, 2) >= 0.03)),

		ERC_tax_local_a1_high1 = cumany(ERC_tax_local_a1 >= ERC_tax_local_a1[year == 1] + 0.05),
		ERC_tax_local_a1_high2 = cumany(ERC_tax_local_a1 >= 0.12),
		ERC_tax_local_a1_hike  = cumany(na2zero(ERC_tax_local_a1 - lag(ERC_tax_local_a1, 2) >= 0.03)),
		
		ERC_tax_GDPState_a1_high1 = cumany(ERC_tax_GDPState_a1 >= ERC_tax_GDPState_a1[year == 1] + 0.05),
		ERC_tax_GDPState_a1_high2 = cumany(ERC_tax_GDPState_a1 >= 0.12),
		ERC_tax_GDPState_a1_hike  = cumany(na2zero(ERC_tax_GDPState_a1 - lag(ERC_tax_GDPState_a1, 2) >= 0.03))
		
		     
		# ERC_tax_PITState_a1_high1  = cumany(ERC_tax_PITState_a1   >= 0.1 ),
		# ERC_tax_PITState_a1_high2  = cumany(ERC_tax_PITState_a1   >= 0.15),
		# ERC_tax_PITState_a1_hike   = cumany(na2zero(ERC_tax_PITState_a1 - lag(ERC_tax_PITState_a1, 2) >= 0.03)),
		# 
		# ERC_tax_salesState_a1_high1 = cumany(ERC_tax_salesState_a1>= 0.1),
		# ERC_tax_salesState_a1_high2 = cumany(ERC_tax_salesState_a1>= 0.15),
		# ERC_tax_salesState_a1_hike  = cumany(na2zero(ERC_tax_salesState_a1 - lag(ERC_tax_salesState_a1, 2) >= 0.03)),
		# 
		# ERC_tax_localState_a1_high1 = cumany(ERC_tax_local_a1 >= 0.1),
		# ERC_tax_localState_a1_high2 = cumany(ERC_tax_local_a1 >= 0.15),
		# ERC_tax_localState_a1_hike  = cumany(na2zero(ERC_tax_local_a1 - lag(ERC_tax_local_a1, 2) >= 0.03))
		
				 ) %>%
	group_by(runname, year) %>%
	dplyr::summarize(
		ERC_tax_PITState_a1_high2 =  sum(ERC_tax_PITState_a1_high2, na.rm = T)/n(),
		ERC_tax_PITState_a1_hike  =  sum(ERC_tax_PITState_a1_hike, na.rm = T)/n(),
		ERC_tax_PITState_a1_high1 =  sum(ERC_tax_PITState_a1_high1, na.rm = T)/n(),
		
		
		ERC_tax_salesState_a1_high2 =  sum(ERC_tax_salesState_a1_high2, na.rm = T)/n(),
		ERC_tax_salesState_a1_hike =  sum(ERC_tax_salesState_a1_hike, na.rm = T)/n(),
		ERC_tax_salesState_a1_high1 =  sum(ERC_tax_salesState_a1_high1, na.rm = T)/n(),
		
		ERC_tax_local_a1_high2   =  sum(ERC_tax_local_a1_high2, na.rm = T)/n(),
		ERC_tax_local_a1_hike    =  sum(ERC_tax_local_a1_hike, na.rm = T)/n(),
		ERC_tax_local_a1_high1   =  sum(ERC_tax_local_a1_high1, na.rm = T)/n(),
		
		ERC_tax_GDPState_a1_high2   =  sum(ERC_tax_GDPState_a1_high2, na.rm = T)/n(),
		ERC_tax_GDPState_a1_hike    =  sum(ERC_tax_GDPState_a1_hike, na.rm = T)/n(),
		ERC_tax_GDPState_a1_high1   =  sum(ERC_tax_GDPState_a1_high1, na.rm = T)/n(),
		
		ERC_tax_PITState_a1.q10 = quantile(ERC_tax_PITState_a1, 0.1,na.rm = T),
		ERC_tax_PITState_a1.q25 = quantile(ERC_tax_PITState_a1, 0.25,na.rm = T),
		ERC_tax_PITState_a1.q50 = quantile(ERC_tax_PITState_a1, 0.5,na.rm = T),
		ERC_tax_PITState_a1.q75 = quantile(ERC_tax_PITState_a1, 0.75,na.rm = T),
		ERC_tax_PITState_a1.q90 = quantile(ERC_tax_PITState_a1, 0.9,na.rm = T),
		
		ERC_tax_salesState_a1.q10 = quantile(ERC_tax_salesState_a1, 0.1,na.rm = T),
		ERC_tax_salesState_a1.q25 = quantile(ERC_tax_salesState_a1, 0.25,na.rm = T),
		ERC_tax_salesState_a1.q50 = quantile(ERC_tax_salesState_a1, 0.5,na.rm = T),
		ERC_tax_salesState_a1.q75 = quantile(ERC_tax_salesState_a1, 0.75,na.rm = T),
		ERC_tax_salesState_a1.q90 = quantile(ERC_tax_salesState_a1, 0.9,na.rm = T),
		
		ERC_tax_local_a1.q10 = quantile(ERC_tax_local_a1, 0.1,na.rm = T),
		ERC_tax_local_a1.q25 = quantile(ERC_tax_local_a1, 0.25,na.rm = T),
		ERC_tax_local_a1.q50 = quantile(ERC_tax_local_a1, 0.5,na.rm = T),
		ERC_tax_local_a1.q75 = quantile(ERC_tax_local_a1, 0.75,na.rm = T),
		ERC_tax_local_a1.q90 = quantile(ERC_tax_local_a1, 0.9,na.rm = T),
		
		ERC_tax_GDPState_a1.q10 = quantile(ERC_tax_GDPState_a1, 0.1,na.rm = T),
		ERC_tax_GDPState_a1.q25 = quantile(ERC_tax_GDPState_a1, 0.25,na.rm = T),
		ERC_tax_GDPState_a1.q50 = quantile(ERC_tax_GDPState_a1, 0.5,na.rm = T),
		ERC_tax_GDPState_a1.q75 = quantile(ERC_tax_GDPState_a1, 0.75,na.rm = T),
		ERC_tax_GDPState_a1.q90 = quantile(ERC_tax_GDPState_a1, 0.9,na.rm = T)
		
		) %>%
	ungroup()





#**********************************************************************
#                         Analysis: Fiscal                                 ####
#**********************************************************************

# Dimensions to compare
  # 1. Types of Stylized governments
  # 2. MS returns vs normal returns
  # 3. Amortization method
  # 4. discout rate

df_risk %<>% select(runname, year, ERC_tax_local_a1_high1, ERC_tax_GDPState_a1_high1, ERC_tax_salesState_a1_high1, ERC_tax_PITState_a1_high1,
									                 ERC_tax_local_a1_high2, ERC_tax_GDPState_a1_high2, ERC_tax_salesState_a1_high2, ERC_tax_PITState_a1_high2,
										               ERC_tax_local_a1_hike,  ERC_tax_GDPState_a1_hike,  ERC_tax_salesState_a1_hike, ERC_tax_PITState_a1_hike) 


## 1. Risk implied by Normally distributed returns and Markov switching returns

df_risk %>% 
	filter(year == 30) %>% 
	select(runname, year, contains("_high1") ) %>% 
	filter(str_detect(runname, "A_"))


df_risk %>% 
	filter(year == 30) %>% 
	select(runname, year, contains("_high2") ) %>% 
	filter(str_detect(runname, "A_"))
	


df_risk %>% 
	filter(year == 30) %>% 
	select(runname, year, contains("_hike") ) %>% 
	filter(str_detect(runname, "A_"))

 # Results: 
  # Risk measures are significantly higher in simulations with synergy between investment returns and economic conditions


# 2. Types of stylized governments

df_risk %>% 
	filter(year == 30, str_detect(runname, "port")) %>% 
	select(runname, year, contains("_high") )

df_risk %>% 
	filter(year == 30, str_detect(runname, "port")) %>% 
	select(runname, year, contains("_hike") )

 # As expected, PIT dominant state faces much greater risk than sales tax state, and benchmark state,
 # especially under very backloaded amortization method

# 3. Amortization method and discount rate

df_risk %>% 
	filter(year == 10, str_detect(runname, "port")) %>% 
	select(runname, year, contains("_high") )

df_risk %>% 
	filter(year == 10, str_detect(runname, "port")) %>% 
	select(runname, year, contains("_hike") )

  # Backloaded amortization method leads to lower contribution risks for government, at the expense of higher risk of severe underfunding,
  # Lower discount rate leads to lower contribution risks for government, at the expense of much higher near-term contributions. 



# Figures 


df_risk %<>%  mutate(
	DiscRate = ifelse(str_detect(runname,"A_"), "Discount rate = 7.5%", "Discount rate = 6.0%"),
	DiscRate = factor(DiscRate, levels = c("Discount rate = 7.5%", "Discount rate = 6.0%")),
	Amort    = ifelse(str_detect(runname,"O30p"), "Slow repayment of UAAL",
										ifelse(str_detect(runname,"O15d"), "Fast repayment of UAAL 15y", "Fast repayment of UAAL 10y")),
	Amort    = factor(Amort, levels = c("Fast repayment of UAAL 10y", "Fast repayment of UAAL 15y", "Slow repayment of UAAL")),
	ReturnDist = ifelse(str_detect(runname,"normal"), "normal", "port70_30"),
	ReturnDist = factor(ReturnDist, levels = c("port70_30", "normal"))
	
	)



## Fig: High ERC as a % of total tax revenue: above 12 percent of payroll
# 
# df_fig1 <- 
# df_risk %>% filter(runname %in% c("A_O10d_port70_30",  "A_O15d_port70_30", "A_O30p_port70_30")) %>% 
#           	select(year, DiscRate, Amort,  year, ERC_tax_PITState_a1_high2, ERC_tax_salesState_a1_high2, ERC_tax_local_a1_high2) 
# 
# df_fig2 <- 
#   df_risk %>% filter(runname %in% c("A_O30pA5_normal", "A_C15d_normal",
# 																	  "B_O30pA5_normal", "B_C15d_normal")) %>% 
# 	select(year, DiscRate, Amort,  year, ERC_tax_local_a1_high2) %>% 
# 	rename(ERC_tax_local_a1_high2_normal = ERC_tax_local_a1_high2)
# df_fig2
# 	
# df_fig <- left_join(df_fig1, df_fig2)
# df_fig
# 
# fig.title <- "Probability of employer contribution as a percentage of total tax revenue \nabove 12 percent at any time up to a given year"
# fig.subtitle <- NULL #"Employer contribution is 5% of total tax revenue in year 1" 
# fig_ERC_tax_high2 <- 
# 	df_fig %>% 
# 	select(year, DiscRate, Amort,  year, ERC_tax_PITState_a1_high2, ERC_tax_salesState_a1_high2, 
# 				                                 ERC_tax_local_a1_high2, ERC_tax_local_a1_high2_normal) %>% 
# 	gather(type, value, -DiscRate, -Amort, -year) %>% 
# 	mutate(type = factor(type, levels = c("ERC_tax_PITState_a1_high2", "ERC_tax_salesState_a1_high2", "ERC_tax_local_a1_high2", "ERC_tax_local_a1_high2_normal"),
# 				                     labels = c("Income tax dominant state \nsimulated returns", "Sales tax dominant state \nsimulated returns", 
# 				                     					  "Baseline:\nconstant tax revenue growth \nsimulated returns", "Baseline:\nconstant tax revenue growth \nnormal returns"))) %>% 
# 	ggplot(aes(x = year, y = 100 * value,
# 						 color = type,
# 						 shape = type)) + 
#   theme_bw() + 
# 	facet_grid(DiscRate ~ Amort) +
# 	geom_line() + 
# 	geom_point(size = 2) + 
# 	geom_hline(yintercept = 100, linetype = 2, size = 1) +
# 	coord_cartesian(ylim = c(0,40)) + 
# 	scale_x_continuous(breaks = c(seq(0, 30, 5))) + 
# 	scale_y_continuous(breaks = seq(0, 500, 5)) + 
# 	scale_color_manual(values = c(color_PIT, color_salesgen, "gray50", "black"),  name = NULL) + 
# 	scale_shape_manual(values = c(15, 16, 17, 18),  name = NULL) +
# 	labs(title = fig.title,
# 			 subtitle = fig.subtitle,
# 			 x = "Year", y = "Probability (%)") + 
# 	theme(axis.text.x = element_text(size = 8)) + 
# 	RIG.theme() + 
# 	theme(legend.position = "bottom")
# fig_ERC_tax_high2
# fig_ERC_tax_high2$data


##  Fig: High ERC as a % of total tax revenue: rising 5 percentage points

df_fig1 <- 
	df_risk %>% filter(runname %in% c("A_O10d_port70_30",  "A_O15d_port70_30", "A_O30p_port70_30")) %>% 
	select(year, DiscRate, Amort,  year, ERC_tax_PITState_a1_high1, ERC_tax_salesState_a1_high1, ERC_tax_local_a1_high1, ERC_tax_GDPState_a1_high1) 

df_fig2 <- 
	df_risk %>% filter(runname %in% c("A_O10d_normal",  "A_O15d_normal", "A_O30p_normal")) %>% 
	select(year, DiscRate, Amort,  year, ERC_tax_local_a1_high1, ERC_tax_GDPState_a1_high1) %>% 
	rename(ERC_tax_local_a1_high1_normal    = ERC_tax_local_a1_high1,
				 ERC_tax_GDPState_a1_high1_normal = ERC_tax_GDPState_a1_high1)
df_fig2

df_fig <- left_join(df_fig1, df_fig2)
df_fig

fig.title <- "Probability of employer contribution as a percentage of total tax revenue \nbeing more than 5 percentage points above the year-1 level at any time up to a given year"
fig.subtitle <- NULL #"Employer contribution is 5% of total tax revenue in year 1" 
fig_ERC_tax_high1 <- 
	df_fig %>% 
	select(year, DiscRate, Amort,  year, 
				 ERC_tax_PITState_a1_high1,
				 ERC_tax_salesState_a1_high1, 
				 ERC_tax_local_a1_high1, 
				 ERC_tax_local_a1_high1_normal,
				 ERC_tax_GDPState_a1_high1, 
				 ERC_tax_GDPState_a1_high1_normal) %>% 
	gather(type, value, -DiscRate, -Amort, -year) %>% 
	mutate(type = factor(type, 
											 levels = c("ERC_tax_local_a1_high1_normal", "ERC_tax_local_a1_high1", "ERC_tax_GDPState_a1_high1_normal", "ERC_tax_GDPState_a1_high1",
											 					  "ERC_tax_salesState_a1_high1",   "ERC_tax_PITState_a1_high1"),
											 labels = c("constant tax revenue growth \nnormal returns",
											 	          "constant tax revenue growth \nsimulated returns", 
											 	          "cyclical tax revenue growth \nnormal returns",
											 	          "cyclical tax revenue growth \nsimulated returns", 
											 	          "Income tax dominant state \nsimulated returns", 
											 					  "Sales tax dominant state \nsimulated returns" 
											 					  )
									
											)) %>% 
	ggplot(aes(x = year, y = 100 * value,
						 color = type,
						 shape = type)) + 
	theme_bw() + 
	facet_grid(DiscRate ~ Amort) +
	geom_line() + 
	geom_point(size = 2) + 
	geom_hline(yintercept = 100, linetype = 2, size = 1) +
	coord_cartesian(ylim = c(0,30)) + 
	scale_x_continuous(breaks = c(seq(0, 30, 5))) + 
	scale_y_continuous(breaks = seq(0, 500, 5)) + 
	scale_color_manual(values = c("grey80", "gray60", "gray40", "black", color_salesgen, color_PIT),  name = NULL) + 
	scale_shape_manual(values = c(15, 16, 17, 18, 19, 20),  name = NULL) +
	labs(title = fig.title,
			 subtitle = fig.subtitle,
			 x = "Year", y = "Probability (%)") + 
	theme(axis.text.x = element_text(size = 8)) + 
	RIG.theme() + 
	theme(legend.position = "bottom")
fig_ERC_tax_high1
fig_ERC_tax_high1$data




## Fig ERC increasing fast relatively to tax revenue 

df_fig1 <- 
	df_risk %>% filter(runname %in% c("A_O10d_port70_30",  "A_O15d_port70_30", "A_O30p_port70_30")) %>% 
	select(year, DiscRate, Amort,  year, ERC_tax_PITState_a1_hike, ERC_tax_salesState_a1_hike, ERC_tax_local_a1_hike, ERC_tax_GDPState_a1_hike) 

df_fig2 <- 
	df_risk %>% filter(runname %in% c("A_O10d_normal",  "A_O15d_normal", "A_O30p_normal")) %>% 
	select(year, DiscRate, Amort,  year, ERC_tax_local_a1_hike, ERC_tax_GDPState_a1_hike) %>% 
	rename(ERC_tax_local_a1_hike_normal    = ERC_tax_local_a1_hike,
				 ERC_tax_GDPState_a1_hike_normal = ERC_tax_GDPState_a1_hike)
df_fig2

df_fig <- left_join(df_fig1, df_fig2)
df_fig


fig.title <- "Probability of employer contribution rising more than 3 percent of total tax revenue\n in a 2-year period at any time up to a given year"
fig.subtitle <- NULL # "Employer contribution is 5% of total tax revenue in year 1" 
fig_ERC_tax_hike <- 
	df_fig %>% 
	select(year, DiscRate, Amort,  year, 
				 ERC_tax_PITState_a1_hike,
				 ERC_tax_salesState_a1_hike, 
				 ERC_tax_local_a1_hike, 
				 ERC_tax_local_a1_hike_normal,
				 ERC_tax_GDPState_a1_hike, 
				 ERC_tax_GDPState_a1_hike_normal) %>% 
	gather(type, value, -DiscRate, -Amort, -year) %>% 
	mutate(type = factor(type, 
											 levels = c("ERC_tax_local_a1_hike_normal", "ERC_tax_local_a1_hike", "ERC_tax_GDPState_a1_hike_normal", "ERC_tax_GDPState_a1_hike",
											 					 "ERC_tax_salesState_a1_hike",   "ERC_tax_PITState_a1_hike"),
											 labels = c("constant tax revenue growth \nnormal returns",
											 					 "constant tax revenue growth \nsimulated returns", 
											 					 "cyclical tax revenue growth \nnormal returns",
											 					 "cyclical tax revenue growth \nsimulated returns", 
											 					 "Income tax dominant state \nsimulated returns", 
											 					 "Sales tax dominant state \nsimulated returns" 
											 )
											 
	)) %>% 
	ggplot(aes(x = year, y = 100 * value,
						 color = type,
						 shape = type)) + 
	theme_bw() + 
	facet_grid(DiscRate ~ Amort) +
	geom_line() + 
	geom_point(size = 2) + 
	geom_hline(yintercept = 100, linetype = 2, size = 1) +
	coord_cartesian(ylim = c(0,60)) + 
	scale_x_continuous(breaks = c(seq(0, 30, 5))) + 
	scale_y_continuous(breaks = seq(0, 500, 10)) + 
	scale_color_manual(values = c("grey80", "gray60", "gray40", "black", color_salesgen, color_PIT),  name = NULL) + 
	scale_shape_manual(values = c(15, 16, 17, 18, 19, 20),  name = NULL) +
	labs(title = fig.title,
			 subtitle = fig.subtitle,
			 x = "Year", y = "Probability (%)") + 
	theme(axis.text.x = element_text(size = 8)) + 
	RIG.theme() + 
	theme(legend.position = "bottom")
fig_ERC_tax_hike
fig_ERC_tax_hike$data


# Saving results

ggsave(paste0(dir_fig_out, "fig_RiskGov_ERC_tax_high1.png"), fig_ERC_tax_high1, width = 8, height = 9 )
#ggsave(paste0(dir_fig_out, "fig_RiskGov_ERC_tax_high2.png"),  fig_ERC_tax_high2,  width = 8, height = 9 )
ggsave(paste0(dir_fig_out, "fig_RiskGov_ERC_tax_hike.png"),   fig_ERC_tax_hike,   width = 8, height = 9 )


Table_ERC_tax_high1 <- 
	df_risk %>% 
	filter(year == 30) %>% 
	select(ReturnDist, DiscRate, Amort, year, contains("_high1") ) %>% 
	arrange(DiscRate, ReturnDist, Amort)
Table_ERC_tax_high1


Table_ERC_tax_high2 <- 
df_risk %>% 
	filter(year == 30) %>% 
	select(ReturnDist, DiscRate, Amort, contains("_high2")) %>% 
	arrange(DiscRate, ReturnDist, Amort)
Table_ERC_tax_high2


Table_ERC_tax_hike <- 
df_risk %>% 
	filter(year == 30) %>% 
	select(ReturnDist, DiscRate, Amort, year, contains("_hike") ) %>% 
	arrange(DiscRate, ReturnDist, Amort)


write.xlsx2(Table_ERC_tax_high1,  file = paste0(dir_fig_out, "Table_RiskGov_ERC_tax_risk.xlsx"), sheetName = "ERC_tax_high1")
write.xlsx2(Table_ERC_tax_high2,  file = paste0(dir_fig_out, "Table_RiskGov_ERC_tax_risk.xlsx"), sheetName = "ERC_tax_high2", append = TRUE)
write.xlsx2(Table_ERC_tax_hike,   file = paste0(dir_fig_out, "Table_RiskGov_ERC_tax_risk.xlsx"), sheetName = "ERC_tax_hike" , append = TRUE)




#**********************************************************************
#                         Analysis: pension finance                ####
#**********************************************************************

df_all.stch <-
	results_all %>%
	filter(sim >= 0, year <= 30)

df_7p5 <- 
results_all %>% 
	filter(runname == "A_O30p_port70_30", sim == 1) %>% 
	select(year, AL_7p5 = AL)


df_all.stch %<>%
	select(runname, sim, year, AL, MA, EEC, PR, ERC_PR) %>%
	left_join(df_7p5) %>% 
	group_by(runname, sim) %>%
	mutate(FR_MA     = 100 * MA / AL_7p5,
				 FR40less   = cumany(FR_MA <= 40),
				 FR100more  = cumany(FR_MA >= 100),
				 FR100more2 = FR_MA >= 100,
				 ERC_high   = cumany(ERC_PR >= 50),
				 ERC_hike   = cumany(na2zero(ERC_PR - lag(ERC_PR, 5) >= 10))) %>%
	group_by(runname, year) %>%
	summarize(FR40less = 100 * sum(FR40less, na.rm = T)/n(),
						FR100more = 100 * sum(FR100more, na.rm = T)/n(),
						FR100more2= 100 * sum(FR100more2, na.rm = T)/n(),
						ERC_high = 100 * sum(ERC_high, na.rm = T)/n(),
						ERC_hike = 100 * sum(ERC_hike, na.rm = T)/n(),
						
						FR.q10   = quantile(FR_MA, 0.1,na.rm = T),
						FR.q25   = quantile(FR_MA, 0.25, na.rm = T),
						FR.q50   = quantile(FR_MA, 0.5, na.rm = T),
						FR.q75   = quantile(FR_MA, 0.75, na.rm = T),
						FR.q90   = quantile(FR_MA, 0.9, na.rm = T),
						
						ERC_PR.q10 = quantile(ERC_PR, 0.1, na.rm = T),
						ERC_PR.q25 = quantile(ERC_PR, 0.25, na.rm = T),
						ERC_PR.q50 = quantile(ERC_PR, 0.5, na.rm = T),
						ERC_PR.q75 = quantile(ERC_PR, 0.75, na.rm = T),
						ERC_PR.q90 = quantile(ERC_PR, 0.9, na.rm = T)
						
						
	) %>%
	ungroup()


df_all.stch %>% filter(year %in% c(1, 10, 20, 30), runname == "A_O30pA5_port70_30")
df_all.stch %>% filter(year %in% c(1, 10, 20, 30), runname == "A_C15d_port70_30" )
df_all.stch %>% filter(year %in% c(1, 10, 20, 30), runname == "A_O30pA5_normal")
df_all.stch %>% filter(year %in% c(1, 10, 20, 30), runname == "A_C15d_normal" )


df_all.stch %>% filter(year %in% c(1, 10, 20, 30), runname == "B_O30pA5_port70_30")
df_all.stch %>% filter(year %in% c(1, 10, 20, 30), runname == "B_C15d_port70_30" )
df_all.stch %>% filter(year %in% c(1, 10, 20, 30), runname == "B_O30pA5_normal")
df_all.stch %>% filter(year %in% c(1, 10, 20, 30), runname == "B_C15d_normal" )



# Payroll growth

results_all %>% 
	filter(sim == 1) %>% 
	select(sim, year, PR, PR.growth)






#*********************************************************************************************************
#  figures 
#*********************************************************************************************************


# # Distribution of funded ratio 
# fig.title <- "Distribution of funded ratios across simulations"
# fig.subtitle <- "Simulated investment returns of a 70/30 portfolio" 
# fig_FRdist <- 
# 	df_all.stch %>% filter(runname %in% c("A_O30pA5_port70_30", "A_C15d_port70_30",
# 																				"B_O30pA5_port70_30", "B_C15d_port70_30"
# 	)) %>% 
# 	select(runname, year, FR.q25, FR.q50, FR.q75) %>% 
# 	gather(type, value, -runname, -year) %>% 
# 	mutate(DiscRate = ifelse(str_detect(runname,"A_"), "Discount rate = 7.5%", "Discount rate = 6.0%"),
# 				 DiscRate = factor(DiscRate, levels = c("Discount rate = 7.5%", "Discount rate = 6.0%")),
# 				 Amort    = ifelse(str_detect(runname,"O30pA5"), "Slow repayment of UAAL", "Fast repayment of UAAL")
# 	) %>% 
# 	mutate(runname = factor(runname, levels = c("A_O30pA5_port70_30", "A_C15d_port70_30"))) %>% 
# 	ggplot(aes(x = year, y = value,
# 						 color = factor(type, levels = c("FR.q75", "FR.q50", "FR.q25")),
# 						 shape = factor(type, levels = c("FR.q75", "FR.q50", "FR.q25"))
# 	)) + theme_bw() + 
# 	facet_grid(DiscRate ~ Amort) +
# 	geom_line() + 
# 	geom_point(size = 2) + 
# 	geom_hline(yintercept = 100, linetype = 2, size = 1) +
# 	coord_cartesian(ylim = c(40,180)) + 
# 	scale_x_continuous(breaks = c(seq(0, 30, 5))) + 
# 	scale_y_continuous(breaks = seq(0, 500, 20)) + 
# 	scale_color_manual(values = c(RIG.green, RIG.blue, RIG.red, "black"), 
# 										 label  = c("75th percentile", "50th percentile", "25th percentile")) + 
# 	scale_shape_manual(values = c(15, 16, 17, 18), 
# 										 label  = c("75th percentile", "50th percentile", "25th percentile")) +
# 	labs(title = fig.title,
# 			 subtitle = fig.subtitle,
# 			 x = "Year", y = "Percent",
# 			 color = NULL, shape = NULL) + 
# 	theme(axis.text.x = element_text(size = 8)) + 
# 	guides(color = guide_legend(keywidth = 2, keyheight = 1))+
# 	theme(legend.position = "bottom") +
# 	RIG.theme()
# fig_FRdist
# fig_FRdist$data
# 
# df_all.stch %>% filter(runname %in% c("A_O30pA5_port70_30", "A_C15d_port70_30") )
# 
# 
# 
# 
# # Distribution of ERC as % Payroll
# fig.title    <- "Distribution of employer contribution as a percentage of payroll across simulations"
# fig.subtitle <- "Simulated investment returns of a 70/30 portfolio"
# fig_ERCdist <- 
# 	df_all.stch %>% filter(runname %in% c("A_O30pA5_port70_30", "A_C15d_port70_30",
# 																				"B_O30pA5_port70_30", "B_C15d_port70_30"
# 	)) %>% 
# 	select(runname, year, ERC_PR.q25, ERC_PR.q50, ERC_PR.q75) %>% 
# 	gather(type, value, -runname, -year) %>% 
# 	mutate(DiscRate = ifelse(str_detect(runname,"A_"), "Discount rate = 7.5%", "Discount rate = 6.0%"),
# 				 DiscRate = factor(DiscRate, levels = c("Discount rate = 7.5%", "Discount rate = 6.0%")),
# 				 Amort    = ifelse(str_detect(runname,"O30pA5"), "Slow repayment of UAAL", "Fast repayment of UAAL")
# 	) %>% 
# 	mutate(runname = factor(runname, levels = c("A_O30pA5_port70_30", "A_C15d_port70_30"))) %>% 
# 	
# 	# mutate(runname = factor(runname, labels = c(lab_s1, lab_s2))) %>%  
# 	ggplot(aes(x = year, y = value,
# 						 color = factor(type, levels = c("ERC_PR.q75", "ERC_PR.q50", "ERC_PR.q25")),
# 						 shape = factor(type, levels = c("ERC_PR.q75", "ERC_PR.q50", "ERC_PR.q25")))) + 
# 	theme_bw() + 
# 	facet_grid(DiscRate ~ Amort) +
# 	geom_line() + 
# 	geom_point(size = 2) + 
# 	geom_hline(yintercept = 100, linetype = 2, size = 1) +
# 	coord_cartesian(ylim = c(0,50)) + 
# 	scale_x_continuous(breaks = c(seq(0, 30, 5))) + 
# 	scale_y_continuous(breaks = seq(0, 500, 5)) + 
# 	scale_color_manual(values = c(RIG.red, RIG.blue, RIG.green, "black"), 
# 										 label  = c("75th percentile", "50th percentile", "25th percentile")) + 
# 	scale_shape_manual(values = c(17, 16, 15, 18), 
# 										 label  = c("75th percentile", "50th percentile", "25th percentile")) +
# 	labs(title = fig.title,
# 			 subtitle = fig.subtitle,
# 			 x = "Year", y = "Percent of payroll",
# 			 color = NULL, shape = NULL) + 
# 	theme(axis.text.x = element_text(size = 8)) + 
# 	guides(color = guide_legend(keywidth = 2, keyheight = 1))+
# 	theme(legend.position = "bottom") +
# 	RIG.theme()
# fig_ERCdist
# 
# 
# # Risk of low funded ratio
# fig.title <- "Probabilities of funded ratio below 40% in any year up to the given year"
# fig.subtitle <- "Simulated investment returns of a 70/30 portfolio"
# fig_FR40less <- 
# 	
# 	df_all.stch %>% filter(runname %in% c("A_O30pA5_port70_30", "A_C15d_port70_30",
# 																				"B_O30pA5_port70_30", "B_C15d_port70_30"
# 	)) %>% 
# 	select(runname, year, FR40less) %>%
# 	gather(type, value, -runname, -year) %>% 
# 	mutate(DiscRate = ifelse(str_detect(runname,"A_"), "Discount rate = 7.5%", "Discount rate = 6.0%"),
# 				 DiscRate = factor(DiscRate, levels = c("Discount rate = 7.5%", "Discount rate = 6.0%")),
# 				 Amort    = ifelse(str_detect(runname,"O30pA5"), "Slow repayment of UAAL", "Fast repayment of UAAL")
# 	) %>% 
# 	mutate(runname = factor(runname, levels = c("A_O30pA5_port70_30", "A_C15d_port70_30"))) %>%
# 	
# 	# mutate(type = factor(type, levels = c("FR75less", "FR60less", "FR40less"), labels = c("75%","60%", "40%" ))) %>% 
# 	#mutate(FR40less.det = 0) %>% 
# 	#gather(variable, value, -year) %>% 
# 	ggplot(aes(x = year, y = value, color = Amort, shape = Amort)) + 
# 	# color = runname, shape = runname)) + 
# 	theme_bw() + 
# 	facet_grid(.~DiscRate) + 
# 	geom_point(size = 2) + 
# 	geom_line() + 
# 	coord_cartesian(ylim = c(0,35)) + 
# 	scale_y_continuous(breaks = seq(0,200, 5)) +
# 	scale_x_continuous(breaks = c(seq(0, 30, 5))) + 
# 	scale_color_manual(values = c(RIG.blue, RIG.green, RIG.red)) + 
# 	scale_shape_manual(values = c(17,16, 15)) +
# 	labs(title = fig.title,
# 			 subtitle = fig.subtitle,
# 			 x = "Year", y = "Probability (%)",
# 			 color = NULL, shape = NULL) + 
# 	guides(color = guide_legend(keywidth = 3, keyheight = 1))+
# 	theme(legend.position = "bottom") + 
# 	RIG.theme()
# fig_FR40less
# fig_FR40less$data %>% filter(year == 2046)
# 
# 
# # Risk of sharp increase in ERC/PR
# fig.title <- "Probability of employer contribution rising more than 10% of payroll \nin a 5-year period at any time up to the given year"
# fig.subtitle <- "Simulated investment returns of a 70/30 portfolio"
# fig_ERChike <- 
# 	
# 	df_all.stch %>% filter(runname %in% c("A_O30pA5_port70_30", "A_C15d_port70_30",
# 																				"B_O30pA5_port70_30", "B_C15d_port70_30"
# 	)) %>% 
# 	select(runname, year, ERC_hike) %>% 
# 	gather(type, value, -runname, -year) %>% 
# 	mutate(DiscRate = ifelse(str_detect(runname,"A_"), "Discount rate = 7.5%", "Discount rate = 6.0%"),
# 				 DiscRate = factor(DiscRate, levels = c("Discount rate = 7.5%", "Discount rate = 6.0%")),
# 				 Amort    = ifelse(str_detect(runname,"O30pA5"), "Slow repayment of UAAL", "Fast repayment of UAAL")
# 	) %>% 
# 	mutate(runname = factor(runname, levels = c("A_O30pA5_port70_30", "A_C15d_port70_30"))) %>%
# 	
# 	ggplot(aes(x = year, y = value, color = Amort, shape = Amort)) + theme_bw() + 
# 	facet_grid(.~DiscRate) + 
# 	geom_point(size = 2) + geom_line() + 
# 	coord_cartesian(ylim = c(0,100)) + 
# 	scale_y_continuous(breaks = seq(0,200, 10)) +
# 	scale_x_continuous(breaks = c(seq(0, 30, 5))) + 
# 	scale_color_manual(values = c(RIG.blue, RIG.green, RIG.red)) + 
# 	scale_shape_manual(values = c(17,16, 15, 18, 19)) +
# 	labs(title = fig.title,
# 			 subtitle = fig.subtitle,
# 			 x = "Year", y = "Probability (%)",
# 			 color = NULL, shape = NULL) + 
# 	guides(color = guide_legend(keywidth = 3, keyheight = 1))+
# 	theme(legend.position = "bottom") + 
# 	RIG.theme()
# fig_ERChike
# fig_ERChike$data %>% filter(year == 2046)
# 

#*********************************************************************************************************
#  Table  
#*********************************************************************************************************

df_all.stch %>% names





Table_risk_DC75_port <- 
	df_all.stch %>%
	select(runname, year, FR40less, ERC_hike, 
				 FR.q10, FR.q25, FR.q50, FR.q75, FR.q90,    
				 ERC_PR.q10, ERC_PR.q25, ERC_PR.q50, ERC_PR.q75, ERC_PR.q90) %>% 
	filter(year %in% c(1, 30), str_detect(runname, "A_"), str_detect(runname, "port")) %>% 
	arrange(year)
Table_risk_DC75_port

Table_risk_DC60_port <-
	df_all.stch %>%
	select(runname, year, FR40less, ERC_hike, 
				 FR.q10, FR.q25, FR.q50, FR.q75, FR.q90,    
				 ERC_PR.q10, ERC_PR.q25, ERC_PR.q50, ERC_PR.q75, ERC_PR.q90) %>% 
	filter(year %in% c(1, 30), str_detect(runname, "B_"), str_detect(runname, "port")) %>% 
	arrange(year)
Table_risk_DC60_port




results_all %>% 
	filter(sim >0, year %in% 1:15, str_detect(runname, "port70_30")) %>%
	group_by(runname, sim) %>% 
	summarise(ERC = sum(ERC/(1.075)^(year - 1))) %>% 
	group_by(runname) %>% 
	summarise(ERC = median(ERC))

results_all %>% 
	filter(sim >0, year %in% 15:30, str_detect(runname, "port70_30")) %>%
	group_by(runname, sim) %>% 
	summarise(ERC = sum(ERC/(1.075)^(year - 1))) %>% 
	group_by(runname) %>% 
	summarise(ERC = median(ERC))

	



Table_risk_DC75 <- 
	df_all.stch %>%
	select(runname, year, FR40less, ERC_hike, 
				 FR.q10, FR.q25, FR.q50, FR.q75, FR.q90,    
				 ERC_PR.q10, ERC_PR.q25, ERC_PR.q50, ERC_PR.q75, ERC_PR.q90) %>% 
	filter(year %in% c(1, 30), str_detect(runname, "A_")) %>% 
	arrange(year)
Table_risk_DC75

Table_risk_DC60 <-
	df_all.stch %>%
	select(runname, year, FR40less, ERC_hike, 
				 FR.q10, FR.q25, FR.q50, FR.q75, FR.q90,    
				 ERC_PR.q10, ERC_PR.q25, ERC_PR.q50, ERC_PR.q75, ERC_PR.q90) %>% 
	filter(year %in% c(1, 30), str_detect(runname, "B_")) %>% 
	arrange(year)
Table_risk_DC60



#*********************************************************************************************************
#  save
#*********************************************************************************************************

# dir_out <- "policyBrief_out/"

write.xlsx2(Table_risk_DC60_port, file = paste0(dir_fig_out, "Table_RiskPen.xlsx"), sheetName = "DC60")
write.xlsx2(Table_risk_DC75_port, file = paste0(dir_fig_out, "Table_RiskPen.xlsx"), sheetName = "DC75", append = TRUE)

# ggsave(fig_FRdist, file  = paste0(dir_fig_out,  "fig_RiskPen_FRdist.png"), width = 8, height = 8 )
# ggsave(fig_ERCdist, file = paste0(dir_fig_out,  "fig_RiskPen_ERCdist.png"), width = 8, height = 8 )
# 
# ggsave(fig_FR40less, file = paste0(dir_fig_out, "fig_RiskPen_FR40less.png"), width = 8, height = 5 )
# ggsave(fig_ERChike,  file = paste0(dir_fig_out, "fig_RiskPen_ERChike.png"),  width = 8, height = 5 )
# 


results_all %>% filter(year == 1, sim ==1) %>% select(runname, FR_MA)