
###########################################################

#### Monte Carlo simulations, analysis of the results #####
####  Definition and production of the output tables  #####

##########################################################


############## Auxiliary packages ###################

library(here)
library(dplyr)
library(data.table)
library(kableExtra)


############## Auxiliary functions ###########

## This function imports a single csv file given as the output of a Monte Carlo
## experiment and do the basic formatting of the data
## Arguments
## x -> a string with the name of the file

import_format <- function(x)
{
  df <- read.csv(here::here('simulation_results',folder,x)) ## import the csv  
  nn <- substr(x,6,nchar(x)-4)  ## remove the word table and extract the simulation number
  df <- cbind(sim_num = rep(nn,nrow(df)),df) # add a column with the simulation number
  setDT(df) ### convert to datatable
  return(df)
}


## This function given a dataset infers what columns are constant, i.e. have the same
## entry for all the rows.
## This is used to store the parameters that are not varying in the MC experiment.
## Arguments
## dt -> an R dataframe

get_constant_columns <- function(df)
{
  setDT(df)
  is_constant <- function(x) length(unique(as.vector(x))) == 1L   
  return(names(df)[sapply(names(df),function(x) is_constant(df[[x]]))])  
}  


########### Select Monte Carlo experiment. Exp_01: BQML finite sample test ###########
folder <- "Exp_01"  ## set the folder
csv_list <- list.files(here::here('simulation_results',folder)) ## read the files list

full_df <- do.call(rbind,lapply(csv_list, import_format))  ## import all the simulation outputs
agg_cols <- colnames(full_df)[17:ncol(full_df)] ## impose aggregation on output variables
group_cols <- colnames(full_df)[c(2:16)] ## group using the parameters of the simulation experiment

## Construct the aggregate dataset containing the MC averages
output <- full_df[,lapply(.SD,function(x) mean(x,na.rm = T)), keyby = group_cols, .SDcols = agg_cols] 


########### Setting up the output tables ################
constant_parameters <- get_constant_columns(output) ## get the parameters that are constant
remove_vars <- c('elapsedTime_BQML','iterations') ## remove the variables that are not relevant (computational performances)
table <-  output[,!(colnames(output) %in% c(constant_parameters,remove_vars)), with =F] ## restyle the output table

parameters_vec_info <- paste0(constant_parameters,': ',output[1,constant_parameters, with = F])
caption_line <- paste(parameters_vec_info, collapse = ', ')


base <- table %>% rename(TSRF = trace_stat_F_BQML) %>% rename(TSRC = trace_stat_C_BQML) %>% rename(MSE = MSE_BQML) 
base[,6:ncol(base)] <- round(base[,6:ncol(base)],4) 


up_table <- base %>% filter(NA_share == 0) 

up_wide <- dcast(up_table, ... ~ idio_cross, value.var = c("TSRC", "TSRF","MSE")) %>% 
           select(-NA_share) %>% 
           select(`T`,n,n_blocks,TSRC_0.1,TSRF_0.1,MSE_0.1,TSRC_0.4,TSRF_0.4,MSE_0.4,TSRC_0.8,TSRF_0.8,MSE_0.8)           
up_wide <- up_wide[order(`T`,n,n_blocks)]
print_up <- up_wide 


minimal_idio_cross <- min(base$idio_cross, na.rm = T)
down_table <-  base %>% filter(idio_cross == minimal_idio_cross, NA_share != 0) 
down_wide <- dcast(down_table, ... ~ NA_share, value.var = c("TSRC", "TSRF","MSE")) %>% 
             select(-idio_cross) %>% 
             select(`T`,n,n_blocks,TSRC_0.5,TSRF_0.5,MSE_0.5,TSRC_0.75,TSRF_0.75,MSE_0.75,TSRC_0.9,TSRF_0.9,MSE_0.9)           
down_wide <- down_wide[order(`T`,n,n_blocks)]
print_down <- down_wide 


options(knitr.table.format = "latex")

out_up <- knitr::kable(print_up,
                       booktabs = T, align = "c", linesep = '',
                       caption = paste0('Missing values share is set to zero. ',caption_line)) %>%
                      kable_styling(latex_options = 'striped') %>%
                      add_header_above(c(" "= 3, "Group 1" = 3, "Group 2" = 3,"Group 3" = 3)) %>%
                      pack_rows(
                        index = c("One global and local factor ($r^C = r^b = 1$)" = 21)
                      )

out_down <- knitr::kable(print_down,
                         booktabs = T, align = "c", linesep = '',
                         caption = paste0('Idiosyncratic cross-correlations is set to ', minimal_idio_cross,'. ', caption_line)) %>%
                          kable_styling(latex_options = 'striped') %>%
                          add_header_above(c(" "= 3, "Group 1" = 3, "Group 2" = 3,"Group 3" = 3)) %>%
                          pack_rows(
                            index = c("One global and local factors ($r^C = r^b = 1$)" = 21)
                          )


writeLines(out_up, here::here('latex_outputs',"Exp_01_1.txt"))
writeLines(out_down, here::here('latex_outputs',"Exp_01_2.txt"))

########### Select Monte Carlo experiment. Exp_02: BE vs BQML with varying block structures ###########

folder <- "Exp_02"
csv_list <- list.files(here::here('simulation_results',folder))

full_df <- do.call(rbind,lapply(csv_list, import_format)) 
agg_cols <- colnames(full_df)[17:ncol(full_df)]
group_cols <- colnames(full_df)[c(2:16)]
output <- full_df[,lapply(.SD,function(x) mean(x,na.rm = T)), keyby = group_cols, .SDcols = agg_cols]


########### Setting up the output tables #######
constant_parameters <- get_constant_columns(full_df) ## get and store the constant parameters
remove_vars <- c('elapsedTime_BQML','elapsedTime_EB') ## remove the variables that are not used

first_table <-  output[,!(colnames(output) %in% c(constant_parameters,remove_vars)), with =F]

parameters_vec_info <- paste0(constant_parameters,': ',output[1,constant_parameters, with = F])
caption_line <- paste(parameters_vec_info, collapse = ', ')

base <- first_table %>% mutate(trace_stat_C_BQML = trace_stat_C_BQML/trace_stat_C_EB) %>%
                        mutate(trace_stat_F_BQML = trace_stat_F_BQML/trace_stat_F_EB) %>%
                        mutate(MSE_BQML = MSE_EB/MSE_BQML) %>%
                        select(-trace_stat_C_EB,-trace_stat_F_EB,-MSE_EB)

tab_homogenous_blocks <- base %>% filter(block_dist == 1) %>% select(-block_dist) %>% rename(TSRF = trace_stat_F_BQML) %>%
                                          rename(TSRC = trace_stat_C_BQML) %>% rename(MSE = MSE_BQML) 
  
tab_heterogeneous_blocks <- base %>% filter(block_dist == 2) %>% select(-block_dist) %>% rename(TSRF = trace_stat_F_BQML) %>%
                                          rename(TSRC = trace_stat_C_BQML) %>% rename(MSE = MSE_BQML)

colnames(tab_homogenous_blocks)[5:7] <- paste0(colnames(tab_homogenous_blocks)[5:7],'_homogenous')
colnames(tab_heterogeneous_blocks)[5:7] <- paste0(colnames(tab_heterogeneous_blocks)[5:7],'_heterogeneous')


print_table <- merge(tab_homogenous_blocks, tab_heterogeneous_blocks,
                     by = c('T','n','r_c','n_blocks'))
print_table <- print_table[order(r_c,`T`,n,n_blocks)] 
print_table[,5:ncol(print_table)] <- round(print_table[,5:ncol(print_table)],4) 
print_table <- print_table %>% select(-r_c)

options(knitr.table.format = "latex")

output_latex <- knitr::kable(print_table,
                             booktabs = T, align = "c", linesep = '',
                             caption = caption_line) %>%
                             kable_styling(latex_options = 'striped') %>%
                             add_header_above(c(" "= 2, "Group 1" = 3, "Group 2" = 3)) %>%
                             pack_rows(
                             index = c("One global factor ($r^C = 1)$" = nrow(print_table)/2,
                                       "Two global factors ($r^C = 2)" = nrow(print_table)/2)
                              )

writeLines(output_latex, here::here('latex_outputs',"Exp_02.txt"))


########### Select Monte Carlo experiment. Exp_03: BE vs BQML with varying cross-correlations and NA share ###########

folder <- "Exp_03" ## set the folder
csv_list <- list.files(here::here('simulation_results',folder)) ## read the files list

full_df <- do.call(rbind,lapply(csv_list, import_format))  ## import all the simulation outputs
agg_cols <- colnames(full_df)[17:ncol(full_df)] ## impose aggregation on output variables
group_cols <- colnames(full_df)[c(2:16)] ## group using the parameters of the simulation experiment

## Construct the aggregate dataset containing the MC averages
output <- full_df[,lapply(.SD,function(x) mean(x,na.rm = T)), keyby = group_cols, .SDcols = agg_cols]


########### Setting up the output tables #######
constant_parameters <- get_constant_columns(output)
remove_vars <- c('elapsedTime_BQML','elapsedTime_EB')
first_table <-  output[,!(colnames(output) %in% c(constant_parameters,remove_vars)), with =F]

parameters_vec_info <- paste0(constant_parameters,': ',output[1,constant_parameters, with = F])
caption_line <- paste(parameters_vec_info, collapse = ', ')


base <- first_table %>% mutate(trace_stat_C_BQML = trace_stat_C_BQML/trace_stat_C_EB) %>%
                        mutate(trace_stat_F_BQML = trace_stat_F_BQML/trace_stat_F_EB) %>%
                        mutate(MSE_BQML = MSE_EB/MSE_BQML) %>%
                        select(-trace_stat_C_EB,-trace_stat_F_EB,-MSE_EB)


a <- base[order(`T`,n,,n_blocks,idio_cross,NA_share)] %>% 
     rename(TSRF = trace_stat_F_BQML) %>%
     rename(TSRC = trace_stat_C_BQML) %>% 
     rename(MSE = MSE_BQML) %>% 
     filter(idio_cross == 0.3) %>%
     dcast(... ~ idio_cross + NA_share, value.var = c('TSRC','TSRF','MSE'))

a <- a[,list(`T`,n,TSRC_0.3_0,TSRF_0.3_0,MSE_0.3_0,TSRC_0.3_0.4,TSRF_0.3_0.4,MSE_0.3_0.6,TSRC_0.3_0.6,TSRF_0.3_0.6,MSE_0.3_0.6)]

b <- base[order(`T`,n,n_blocks,idio_cross,NA_share)] %>% 
     rename(TSRF = trace_stat_F_BQML) %>%
     rename(TSRC = trace_stat_C_BQML) %>% 
     rename(MSE = MSE_BQML) %>% 
     filter(idio_cross == 0.6) %>%
     dcast(... ~ idio_cross + NA_share, value.var = c('TSRC','TSRF','MSE'))

b <- b[,list(`T`,n,TSRC_0.6_0,TSRF_0.6_0,MSE_0.6_0,TSRC_0.6_0.4,TSRF_0.6_0.4,MSE_0.6_0.4,TSRC_0.6_0.6,TSRF_0.6_0.6,MSE_0.6_0.6)]

a[,3:ncol(a)] <- round(a[,3:ncol(a)],4)
b[,3:ncol(b)] <- round(b[,3:ncol(b)],4)

options(knitr.table.format = "latex")

a_out <- knitr::kable(a,
                    booktabs = T, align = "c", linesep = '',
                    caption = paste0('Idiosyncratic cross-correlation is 0.3. ',caption_line)) %>%
                    kable_styling(latex_options = 'striped') %>%
                    add_header_above(c(" "= 2, "tau3 delta0" = 3, "tau3 delta4" = 3,  "tau3 delta6" = 3)) 

b_out <- knitr::kable(b,
                      booktabs = T, align = "c", linesep = '',
                      caption = paste0('Idiosyncratic cross-correlation is 0.6. ',caption_line)) %>%
                      kable_styling(latex_options = 'striped') %>%
                      add_header_above(c(" "= 2, "tau6 delta0" = 3, "tau6 delta4" = 3,  "tau6 delta6" = 3)) 


writeLines(a_out, here::here('latex_outputs',"Exp_03_1.txt"))
writeLines(b_out, here::here('latex_outputs',"Exp_03_2.txt"))


########### Select Monte Carlo experiment. Exp_04: BQML and BE under fat-tailed distributions of the input series #########

folder <- "Exp_04" ## set the folder
csv_list <- list.files(here::here('simulation_results',folder)) ## read the files list

full_df <- do.call(rbind,lapply(csv_list, import_format))  ## import all the simulation outputs
agg_cols <- colnames(full_df)[17:ncol(full_df)] ## impose aggregation on output variables
group_cols <- colnames(full_df)[c(2:16)] ## group using the parameters of the simulation experiment

## Construct the aggregate dataset containing the MC averages
output <- full_df[,lapply(.SD,function(x) mean(x,na.rm = T)), keyby = group_cols, .SDcols = agg_cols]


########### Setting up the output tables #######
constant_parameters <- get_constant_columns(output)
remove_vars <- c('elapsedTime_BQML','elapsedTime_EB')
first_table <-  output[,!(colnames(output) %in% c(constant_parameters,remove_vars)), with =F]

parameters_vec_info <- paste0(constant_parameters,': ',output[1,constant_parameters, with = F])
caption_line <- paste(parameters_vec_info, collapse = ', ')


base_compare <- first_table %>% mutate(trace_stat_C_BQML = trace_stat_C_BQML/trace_stat_C_EB) %>%
                                mutate(trace_stat_F_BQML = trace_stat_F_BQML/trace_stat_F_EB) %>%
                                mutate(MSE_BQML = MSE_EB/MSE_BQML) %>%
                                select(-trace_stat_C_EB,-trace_stat_F_EB,-MSE_EB,-MSE_BQML) %>%
                                rename(TSRF_c = trace_stat_F_BQML) %>% rename(TSRC_c = trace_stat_C_BQML) 


base <- first_table %>% select(-trace_stat_C_EB,-trace_stat_F_EB,-MSE_EB,-MSE_BQML)  %>%
                        rename(TSRF = trace_stat_F_BQML) %>% rename(TSRC = trace_stat_C_BQML)


base <- merge(base,base_compare, by = c('T','n','idio_dist','idio_cross','NA_share'))
base[,6:ncol(base)] <- round(base[,6:ncol(base)],4)


###
a_t1_d0 <- base[order(`T`,n,idio_cross,NA_share,idio_dist)] %>% filter(idio_dist %in% c('laplace','student_3')) %>%
                                                                filter(idio_cross == 0.1,NA_share == 0) %>%
          dcast(... ~ idio_dist, value.var = c('TSRC','TSRF','TSRC_c','TSRF_c'))

a_t1_d0 <- a_t1_d0[,list(`T`,n,TSRC_laplace,TSRF_laplace,TSRC_student_3,TSRF_student_3,
                               TSRC_c_laplace,TSRF_c_laplace,TSRC_c_student_3,TSRF_c_student_3)]


###
a_t1_d5 <- base[order(`T`,n,idio_cross,NA_share,idio_dist)] %>% filter(idio_dist %in% c('laplace','student_3')) %>%
                        filter(idio_cross == 0.1,NA_share == 0.5) %>%
                        dcast(... ~ idio_dist, value.var = c('TSRC','TSRF','TSRC_c','TSRF_c'))

a_t1_d5 <- a_t1_d5[,list(`T`,n,TSRC_laplace,TSRF_laplace,TSRC_student_3,TSRF_student_3,
                         TSRC_c_laplace,TSRF_c_laplace,TSRC_c_student_3,TSRF_c_student_3)]


###
a_t5_d0 <- base[order(`T`,n,idio_cross,NA_share,idio_dist)] %>% filter(idio_dist %in% c('laplace','student_3')) %>%
          filter(idio_cross == 0.5,NA_share == 0) %>%
          dcast(... ~ idio_dist, value.var = c('TSRC','TSRF','TSRC_c','TSRF_c'))

a_t5_d0 <- a_t5_d0[,list(`T`,n,TSRC_laplace,TSRF_laplace,TSRC_student_3,TSRF_student_3,
                         TSRC_c_laplace,TSRF_c_laplace,TSRC_c_student_3,TSRF_c_student_3)]


###
a_t5_d5 <- base[order(`T`,n,idio_cross,NA_share,idio_dist)] %>% filter(idio_dist %in% c('laplace','student_3')) %>%
                filter(idio_cross == 0.5,NA_share == 0.5) %>%
                dcast(... ~ idio_dist, value.var = c('TSRC','TSRF','TSRC_c','TSRF_c'))


a_t5_d5 <- a_t5_d5[,list(`T`,n,TSRC_laplace,TSRF_laplace,TSRC_student_3,TSRF_student_3,
                         TSRC_c_laplace,TSRF_c_laplace,TSRC_c_student_3,TSRF_c_student_3)]


options(knitr.table.format = "latex")

### Preparing output tables in latex format
out_t1_d0 <- knitr::kable(a_t1_d0,
                          booktabs = T, align = "c", linesep = '',
                          caption = paste0('Idiosyncratic cross-correlations is 0.1. No missing values. ',caption_line)) %>%
                          kable_styling(latex_options = 'striped') %>%
                          add_header_above(c(" "= 2, "Laplace" = 2, "Student 3" = 2, "Laplace" = 2, "Student 3" = 2)) %>%
                          pack_rows(
                                    index = c("tau = 0.1 delta = 0" = nrow(a_t1_d0))
                                    )

out_t1_d5 <- knitr::kable(a_t1_d5,
                        booktabs = T, align = "c", linesep = '',
                        caption = paste0('Idiosyncratic cross-correlations is 0.1. Missing values share is 0.5. ', caption_line)) %>%
                        kable_styling(latex_options = 'striped') %>%
                        add_header_above(c(" "= 2, "Laplace" = 2, "Student 3" = 2, "Laplace" = 2, "Student 3" = 2)) %>%
                        pack_rows(
                                  index = c("tau = 0.1 delta = 0.5" = nrow(a_t1_d0))
                                  )

out_t5_d0 <- knitr::kable(a_t5_d0,
                          booktabs = T, align = "c", linesep = '',
                          caption = paste0('Idiosyncratic cross-correlations is 0.5. No missing values. ', caption_line)) %>%
                          kable_styling(latex_options = 'striped') %>%
                          add_header_above(c(" "= 2, "Laplace" = 2, "Student 3" = 2, "Laplace" = 2, "Student 3" = 2)) %>%
                          pack_rows(
                                    index = c("tau = 0.5 delta = 0" = nrow(a_t1_d0))
                                    )

out_t5_d5 <- knitr::kable(a_t5_d5,
                        booktabs = T, align = "c", linesep = '',
                        caption = paste0('Idiosyncratic cross-correlations is 0.5. Missing values share is 0.5. ',caption_line)) %>%
                        kable_styling(latex_options = 'striped') %>%
                        add_header_above(c(" "= 2, "Laplace" = 2, "Student 3" = 2, "Laplace" = 2, "Student 3" = 2)) %>%
                        pack_rows(
                                  index = c("tau = 0.5 delta = 0.5" = nrow(a_t1_d0))
                                  )


writeLines(out_t1_d0, here::here('latex_outputs',"Exp_04_1.txt"))
writeLines(out_t1_d5, here::here('latex_outputs',"Exp_04_2.txt"))
writeLines(out_t5_d0, here::here('latex_outputs',"Exp_04_3.txt"))
writeLines(out_t5_d5, here::here('latex_outputs',"Exp_04_4.txt"))