##### Function to produce the summary statistics and plots from the dataset
##### Inputs: base_filt <- dataset to be analyzed
#####         folder <- folder path where to save the outcomes

#### The statistics are collected in tables saved in latex format. Ready to be copied
#### and inserted in a latex main.

if (!require("pacman")) install.packages("pacman")


pacman::p_load(data.table)
pacman::p_load(tidyr)
pacman::p_load(dplyr)
pacman::p_load(here)
pacman::p_load(matrixStats)


data_summary_function <- function
(
  base_filt, ## dataset with only firm and country dimension, zeros when missing transaction
  folder = 'data_summary' ## folder where text elements, summary dataset and plots are saved
)
{
  ## Define a function converting the NAs to zeros
  f_zero = function(DT) 
  {
    for (i in names(DT))
      DT[is.na(get(i)), (i):=0]
  }
  ## Define a function converting the zeros to NAs
  f_nan = function(DT) 
  {
    for (i in names(DT))
      DT[get(i)==0, (i):=NA]
  }
  
  ##### Create the folders where to store the outcomes (if not existing) #####
  dir.create(here::here('data_preparation_sampling'), showWarnings = FALSE)
  dir.create(here::here('data_preparation_sampling',folder), showWarnings = FALSE)
  
  ## Preliminary operation: identify and store the colunm names of the export value columns
  export_cols <- colnames(base_filt)[3:ncol(base_filt)]
  
  ## Force numeric format for columns containing export values
  base_filt[,(export_cols) := lapply(.SD, function(x) as.numeric(x)), .SDcols = export_cols]
  
  ################################## Calculate quarterly GRs and volatility #######################################

  f_nan(base_filt) ## substituting NAs with zeros (at this stage they are equivalent)
  mat_base <- as.matrix(base_filt[,3:ncol(base_filt)]) ## extract the numeric matrix
  gr_mat <- log(1+mat_base[,5:ncol(mat_base)]) - log(1+mat_base[,1:(ncol(mat_base)-4)]) ## calculate GR
  gr_base_q <- setDT(as.data.frame(gr_mat))  ## convert to data.table
  
  base_filt_q <- base_filt[rowSums(!is.na(gr_base_q), na.rm = T)>1] ## remove empty obs in GR from original dataset
  gr_base_q <- gr_base_q[rowSums(!is.na(gr_base_q), na.rm = T)>1] ## remove empty obs in GR from GR dataset
  gr_base_q <- cbind(base_filt_q[,1:2, with = F],gr_base_q) ## insert siren pyod indicators
  
  rm(mat_base)
  rm(gr_mat)
  
  
  ## calculate the log growth rates of the firm ##
  siren_agg_q <- base_filt_q[, lapply(.SD, function(x) sum(as.numeric(x), na.rm = T)),
                             by =c('SIREN'), 
                             .SDcols = export_cols]
  
  mat_siren_base <- as.matrix(siren_agg_q[,3:ncol(siren_agg_q)])
  gr_siren <- log(1+mat_siren_base[,5:ncol(mat_siren_base)]) - log(1+mat_siren_base[,1:(ncol(mat_siren_base)-4)])
  gr_siren_q <- cbind(siren_agg_q[,1:2, with = F],gr_siren)
  
  ## Compute the weights for the weighted sum of the average
  weight_siren <-  base_filt_q[ ,
                                paste0(export_cols,'w') := lapply(.SD, function(x) as.numeric(x)/sum(as.numeric(x), na.rm = T)),
                                by = c('SIREN'), .SDcols = export_cols][,c('PYOD','SIREN',paste0(export_cols,'w')), with = F]
  
  ## Weighted logGR i.e. logGR = w_{t-1}(log(y_t)-log(y_{t-1})
  wgt_log_gr <-  gr_base_q[,3:98,with = F]*weight_siren[ ,3:98, with =F]
  
  wgt_log_gr <- cbind(gr_base_q[,1:2, with =F],wgt_log_gr)
  col_q  <- colnames(wgt_log_gr)[3:98]
  wgt_log_gr_siren <- wgt_log_gr[,lapply(.SD, function(x)sum(x, na.rm = T)),  keyby = c('SIREN'),.SDcols = col_q]
  f_nan(wgt_log_gr_siren)
  
  ## Create an object with weighted volatilities (quarterly time unit) per siren

  wgt_vol_q <- cbind(wgt_log_gr_siren[, c(1), with = F],
                     WGT_VOL_Q = rowVars(as.matrix(wgt_log_gr_siren[,-c(1), with = F]), na.rm = T))
  vol_q <- cbind(gr_siren_q[, c(1), with = F], VOL_Q = rowVars(as.matrix(gr_siren_q[,-c(1), with = F]), na.rm = T))
  
  setkey(vol_q,SIREN)
  setkey(wgt_vol_q,SIREN)
  
  vol_compare <- merge(vol_q, wgt_vol_q, all = T)


  ################################## Calculate yearly GRs and volatility #######################################
  
  ############################## Transform data in yearly base #################################
  
  base_filt_y <- base_filt[,1:2]
  colnames(base_filt_y)[1:2] <- colnames(base_filt)[1:2]
  n_years <- floor((ncol(base_filt)-2)/4)
  for(s in 1:n_years)
  {
    base_filt_y <- cbind(base_filt_y,rowSums(base_filt[,2+(4*(s-1)+1):(4*(s-1)+4),with = F], na.rm = T))
    colnames(base_filt_y)[2+s] <- substr(colnames(base_filt)[2+(4*(s-1)+1)],1,4)
  } 

  f_nan(base_filt_y) ## substituting NAs with zeros (at this stage they are equivalent)
  
  mat_base <- as.matrix(base_filt_y[,3:ncol(base_filt_y)]) ## extract the numeric matrix
  gr_mat <- log(1+mat_base[,2:ncol(mat_base)]) - log(1+mat_base[,1:(ncol(mat_base)-1)]) ## calculate GR
  gr_base <- setDT(as.data.frame(gr_mat))  ## convert to data.table
  
  base_filt_y <- base_filt_y[rowSums(!is.na(gr_base), na.rm = T)>1] ## remove empty obs in GR from original dataset
  gr_base <- gr_base[rowSums(!is.na(gr_base), na.rm = T)>1] ## remove empty obs in GR from GR dataset
  gr_base <- cbind(base_filt_y[,1:2, with = F],gr_base) ## insert siren pyod indicators
  
  
  ## Compute the log growth rate of the firm (yearly)
  siren_agg_y <- base_filt_y[, lapply(.SD, function(x) sum(as.numeric(x), na.rm = T)),
                             by =c('SIREN'), 
                             .SDcols = as.character(1993:2017)]
  
  mat_siren_base <- as.matrix(siren_agg_y[,3:ncol(siren_agg_y)])
  gr_siren <- log(1+mat_siren_base[,2:ncol(mat_siren_base)]) - log(1+mat_siren_base[,1:(ncol(mat_siren_base)-1)])
  gr_siren_y <- cbind(siren_agg_y[,1:2, with = F],gr_siren)
  
  
  ## Compute the weights for the weighted sum of the average
  cols_y <- colnames(base_filt_y)[-c(1,2)]
  
  weight_siren_y <-  base_filt_y[ ,
                                  paste0(cols_y,'w') :=lapply(.SD, function(x) x/sum(x, na.rm = T)),
                                  by = c('SIREN'), .SDcols = cols_y][,c('PYOD','SIREN',paste0(cols_y,'w')), with = F]
  
  ## Weighted logGR i.e. logGR = w_{t-1}(log(y_t)-log(y_{t-1}) ##
  wgt_log_gr <-  gr_base[,3:26,with = F]*weight_siren_y[ ,3:26, with =F]
  
  wgt_log_gr <- cbind(gr_base[,1:2, with =F],wgt_log_gr)
  col_q  <- colnames(wgt_log_gr)[3:26]
  wgt_log_gr_siren <- wgt_log_gr[,lapply(.SD, function(x)sum(x, na.rm = T)),  keyby = c('SIREN'),.SDcols = col_q]
  f_nan(wgt_log_gr_siren)
  
  ## Create an object with weighted volatilities (quarterly time unit) per siren
  
  wgt_vol_y <- cbind(wgt_log_gr_siren[, c(1), with = F],
                     WGT_VOL_Y = rowVars(as.matrix(wgt_log_gr_siren[,-c(1), with = F]), na.rm = T))
  
  vol_y <- cbind(gr_siren_y[, c(1), with = F], VOL_Y = rowVars(as.matrix(gr_siren_y[,-c(1), with = F]), na.rm = T))
  
  setkey(vol_y,SIREN)
  setkey(wgt_vol_y,SIREN)
  
  vol_compare <- merge(vol_compare,merge(vol_y, wgt_vol_y, all = T), all = T)
  
  rm(gr_mat)
  rm(wgt_log_gr_siren)
  rm(weight_siren_y)
  rm(wgt_log_gr)
  
  ######################################## Data summary ########################################
  
  #################################### Auxiliary functions #####################################
  
  H_in <- function(x) ## given a vector of shares calculate the Herfindal index
  {
    sum(x^2, na.rm = T)
  }

  DIV_in <- function(x) ## given a vector of shares calculate a diversification index as the opposite of the HH index
  {
  ind <- 1/sum(x^2, na.rm = T)  
  out <- ifelse(is.infinite(ind),NA,ind)
  return(out)
  }  
  
  S_in <- function(x) ## given a vector of shares calculate the Shannon index
  {
    x[is.na(x)] <- 0
    y <- -x%*%log(x)
    return(y)
  }
  
  
  ############# Produce yearly statistics per siren  ###########
  
  tab_sum <- cbind(base_filt_y[,1:2],
                   VAL = rowSums(base_filt_y[,-c('PYOD','SIREN')], na.rm = T))
  
  tab_sum[, NUM := .N, keyby = c('SIREN')]   ## number of destinations per firm
  tab_sum[, LOG_NUM := log(NUM)] ## log number of destinations per firm
  tab_sum[, SH := VAL/sum(VAL, na.rm =T), keyby = c('SIREN')] ## share of value of each flow per siren
  tab_sum[, S_TOP :=  max(SH), keyby = c('SIREN')] ## top share
  tab_sum[, S_IN :=  S_in(SH), keyby = c('SIREN')] ## shannon index of the shares
  tab_sum[, H_IN := H_in(SH), by = c('SIREN')] ## herfindal index (concentration) of the shares
  tab_sum[, DIV_IN := DIV_in(SH), by = c('SIREN')] ## diversification index of the shares
  tab_sum <- unique(tab_sum[,-c('VAL','SH','PYOD'),with = F])  ## compute the firm level statistics
  
  setkey(tab_sum,SIREN)
  setkey(wgt_vol_y,SIREN)
  setkey(wgt_vol_q,SIREN)
  
  tab_sum <- merge(tab_sum,wgt_vol_y, all = T)
  tab_sum <- merge(tab_sum,wgt_vol_q, all = T)
  rm(wgt_vol_q)
  rm(wgt_vol_y)
  
  ############ Produce some data summary - Summary year by year ##############
  ################################### Firm level statisics ####################################
  tab_sum_calc <- base_filt_y[,1:(n_years+2)]
  tab_sum_y <- base_filt_y[,1:2]

  cols_y <- colnames(tab_sum_calc)[3:ncol(tab_sum_calc)]
  
  ### Mean and median number of destinations in the whole time span - yearly frewquency
  num_sum_y <- tab_sum_calc[, paste0('NUM_',cols_y) := lapply(.SD,function(x) sum(!is.na(x))), keyby = c('SIREN'), .SDcols = cols_y]
  tab_sum_y <- cbind(tab_sum_y, MEAN_NUM = rowMeans(num_sum_y[,paste0('NUM_',cols_y), with = F],na.rm = T))
  tab_sum_y[, LOG_MEAN_NUM := log(MEAN_NUM)]
  tab_sum_y <- cbind(tab_sum_y, MEDIAN_NUM = rowMedians(as.matrix(num_sum_y[,paste0('NUM_',cols_y), with = F]),na.rm = T))
  tab_sum_y[, LOG_MED_NUM := log(MEDIAN_NUM)]
  
  ### AVG size of the firm export value (in the whole time span - yearly frequency)
  siren_size <- cbind(tab_sum_calc[,1:2],AVG_SIZE = rowMeans(tab_sum_calc[,3:(n_years+2)], na.rm = T))
  siren_size <- unique(siren_size[,AVG_SIZE := mean(AVG_SIZE), keyby = c('SIREN')][,-c('PYOD')])
  
  ### Calculating yearly diversifiaction and concentration indices
  tab_sum_calc <- tab_sum_calc[,1:(n_years+2)]
  
  ## Calculating the shares
  sh_sum_y <- tab_sum_calc[,paste0('SH_',cols_y) := lapply(.SD,function(x) x/sum(x, na.rm =T)),
                           keyby = c('SIREN'),
                           .SDcols = cols_y]
  sh_sum_y <- sh_sum_y[,c('SIREN','PYOD',paste0('SH_',cols_y)), with = F]
  
  f_zero(sh_sum_y)
  
  S_in_mod <- function(x){as.vector(S_in(x[x!=0]))*rep(1,length(x))} ## shannon index function
  H_in_mod <- function(x){H_in(x)*rep(1,length(x))} ## herfindahl index function
  DIV_in_mod <- function(x){DIV_in(x)*rep(1,length(x))} ## herfindahl index function
  
  ### Year by year Shannon index  
  sh_sum_y <- sh_sum_y[,c('SIREN','PYOD',paste0('SH_',cols_y)), with = F][,paste0('S_IN',cols_y) := lapply(.SD,function(x)S_in_mod(x)),
                                                                          keyby = c('SIREN'),
                                                                          .SDcols = paste0('SH_',cols_y)]
  
  ## Time span (yearly) mean and median of the Shannon index
  tab_sum_y <- cbind(tab_sum_y,S_IN_MEAN = rowMeans(sh_sum_y[,paste0('S_IN',cols_y),with = F], na.rm = T))
  tab_sum_y <- cbind(tab_sum_y,S_IN_MED = rowMedians(as.matrix(sh_sum_y[,paste0('S_IN',cols_y),with = F]), na.rm = T))
  
 
  ### Year by year HHindex  
  sh_sum_y <- sh_sum_y[,c('SIREN','PYOD',paste0('SH_',cols_y)), with = F][,paste0('H_IN',cols_y) := lapply(.SD,function(x)H_in_mod(x)),
                                                                          keyby = c('SIREN'),
                                                                          .SDcols = paste0('SH_',cols_y)]
  
  
  ## Time span (yearly) mean and median of the Herfindahl index
  tab_sum_y <- cbind(tab_sum_y,H_IN_MEAN = rowMeans(sh_sum_y[,paste0('H_IN',cols_y),with = F], na.rm = T))
  tab_sum_y <- cbind(tab_sum_y,H_IN_MED = rowMedians(as.matrix(sh_sum_y[,paste0('H_IN',cols_y),with = F]), na.rm = T))
  
  
  ### Year by year HHindex  
  sh_sum_y <- sh_sum_y[,c('SIREN','PYOD',paste0('SH_',cols_y)), with = F][,paste0('DIV_IN',cols_y) := lapply(.SD,function(x)DIV_in_mod(x)),
                                                                          keyby = c('SIREN'),
                                                                          .SDcols = paste0('SH_',cols_y)]
  
  
  ## Time span (yearly) mean and median of the Herfindahl index
  tab_sum_y <- cbind(tab_sum_y,DIV_IN_MEAN = rowMeans(sh_sum_y[,paste0('DIV_IN',cols_y),with = F], na.rm = T))
  tab_sum_y <- cbind(tab_sum_y,DIV_IN_MED = rowMedians(as.matrix(sh_sum_y[,paste0('DIV_IN',cols_y),with = F]), na.rm = T))
  

  
  ## Year by year Top Share  
  sh_sum_y <- sh_sum_y[,c('SIREN','PYOD',paste0('SH_',cols_y)), with = F][,paste0('S_TOP',cols_y) := lapply(.SD,function(x)max(x)),
                                                                          keyby = c('SIREN'),
                                                                          .SDcols = paste0('SH_',cols_y)]
  ## Time span (yearly) mean and median of the top share
  tab_sum_y <- cbind(tab_sum_y,S_TOP_MEAN = rowMeans(sh_sum_y[,paste0('S_TOP',cols_y),with = F], na.rm = T))
  tab_sum_y <- cbind(tab_sum_y,S_TOP_MED = rowMedians(as.matrix(sh_sum_y[,paste0('S_TOP',cols_y),with = F]), na.rm = T))
  
  
  ## Cleaning and merging with firm level statistics
  tab_sum_y <- unique(tab_sum_y[,-c('PYOD')])
  
  colnames(tab_sum_y)[-1] <- paste0(colnames(tab_sum_y)[-1],'_Y')
  
  setkey(tab_sum_y,SIREN)
  setkey(tab_sum,SIREN)
  
  tab_sum <- merge(tab_sum,tab_sum_y, all = T)
  
  rm(tab_sum_y)
  
  
  ##########  Summary quarter by quarter ###########
  
  ################################### Summary quarter by quarter ##################################
  
  tab_sum_calc <- base_filt[,c('SIREN','PYOD',export_cols),with  = F]
  tab_sum_q <- base_filt[,1:2]
  
  cols_q <- colnames(tab_sum_calc)[3:ncol(tab_sum_calc)]
  
  ### Mean and median number of destinations in the whole time span - quarterly frewquency
  num_sum_q <- tab_sum_calc[, paste0('NUM_',cols_q) := lapply(.SD,function(x) sum(!is.na(x))), keyby = c('SIREN'), .SDcols = cols_q]
  tab_sum_q <- cbind(tab_sum_q, MEAN_NUM = rowMeans(num_sum_q[,paste0('NUM_',cols_q), with = F],na.rm = T))
  tab_sum_q[, LOG_MEAN_NUM := log(MEAN_NUM)]
  tab_sum_q <- cbind(tab_sum_q, MEDIAN_NUM = rowMedians(as.matrix(num_sum_q[,paste0('NUM_',cols_q), with = F]),na.rm = T))
  tab_sum_q[, LOG_MED_NUM := log(MEDIAN_NUM)]
  

  tab_sum_calc <- tab_sum_calc[,c('SIREN','PYOD',export_cols),with  = F]
  
  ## Calculating the shares
  sh_sum_q <- tab_sum_calc[,paste0('SH_',cols_q) := lapply(.SD,function(x) x/sum(x, na.rm =T)),
                           keyby = c('SIREN'),
                           .SDcols = cols_q]
  sh_sum_q <- sh_sum_q[,c('SIREN','PYOD',paste0('SH_',cols_q)), with = F]
  
  f_zero(sh_sum_q) ## replace NA with zeros
  

  ## Year by year Shannon index   
  sh_sum_q <- sh_sum_q[,c('SIREN','PYOD',paste0('SH_',cols_q)), with = F][,paste0('S_IN',cols_q) := lapply(.SD,function(x)S_in_mod(x)),
                                                                          keyby = c('SIREN'),
                                                                          .SDcols = paste0('SH_',cols_q)]
  
  ## Time span (quarterly) mean and median of the Shannon index
  tab_sum_q <- cbind(tab_sum_q, S_IN_MEAN = rowMeans(sh_sum_q[,paste0('S_IN',cols_q),with = F], na.rm = T))
  tab_sum_q <- cbind(tab_sum_q, S_IN_MED = rowMedians(as.matrix(sh_sum_q[,paste0('S_IN',cols_q),with = F]), na.rm = T))
  
  
  ## Year by year HH index   
  sh_sum_q <- sh_sum_q[,c('SIREN','PYOD',paste0('SH_',cols_q)), with = F][,paste0('H_IN',cols_q) := lapply(.SD,function(x)H_in_mod(x)),
                                                                          keyby = c('SIREN'),
                                                                          .SDcols = paste0('SH_',cols_q)]
  
  
  ## Time span (quarterly) mean and median of the Herfindahl index
  tab_sum_q <- cbind(tab_sum_q,H_IN_MEAN = rowMeans(sh_sum_q[,paste0('H_IN',cols_q),with = F], na.rm = T))
  tab_sum_q <- cbind(tab_sum_q,H_IN_MED = rowMedians(as.matrix(sh_sum_q[,paste0('H_IN',cols_q),with = F]), na.rm = T))
  
  ## Year by year diversification index   
  sh_sum_q <- sh_sum_q[,c('SIREN','PYOD',paste0('SH_',cols_q)), with = F][,paste0('DIV_IN',cols_q) := lapply(.SD,function(x)DIV_in_mod(x)),
                                                                          keyby = c('SIREN'),
                                                                          .SDcols = paste0('SH_',cols_q)]
  
  
  ## Time span (quarterly) mean and median of the Herfindahl index
  tab_sum_q <- cbind(tab_sum_q,DIV_IN_MEAN = rowMeans(sh_sum_q[,paste0('DIV_IN',cols_q),with = F], na.rm = T))
  tab_sum_q <- cbind(tab_sum_q,DIV_IN_MED = rowMedians(as.matrix(sh_sum_q[,paste0('DIV_IN',cols_q),with = F]), na.rm = T))
  

  
  ## Year by year top share   
  sh_sum_q <- sh_sum_q[,c('SIREN','PYOD',paste0('SH_',cols_q)), with = F][,paste0('S_TOP',cols_q) := lapply(.SD,function(x)max(x)),
                                                                          keyby = c('SIREN'),
                                                                          .SDcols = paste0('SH_',cols_q)]
  ## Time span (quarterly) mean and median of the top share
  tab_sum_q <- cbind(tab_sum_q,S_TOP_MEAN = rowMeans(sh_sum_q[,paste0('S_TOP',cols_q),with = F], na.rm = T))
  tab_sum_q <- cbind(tab_sum_q,S_TOP_MED = rowMedians(as.matrix(sh_sum_q[,paste0('S_TOP',cols_q),with = F]), na.rm = T))
  
  
  ## Clean and merge the statistics on the quarterly frequency
  tab_sum_q <- unique(tab_sum_q[,-c('PYOD')])
  
  colnames(tab_sum_q)[-1] <- paste0(colnames(tab_sum_q)[-1],'_Q')
  
  setkey(tab_sum_q,SIREN)
  setkey(tab_sum,SIREN)
  
  tab_sum <- merge(tab_sum,tab_sum_q, all = T)
  
  setkey(siren_size,SIREN)
  tab_sum <- merge(tab_sum,siren_size, all = T)
  tab_sum <- merge(tab_sum, vol_compare[,c('SIREN','VOL_Q','VOL_Y'), with = F])
  
  
  cols <- colnames(tab_sum)[-1]
  cor_mat <- cor(tab_sum[,lapply(.SD, as.numeric),.SDcols = cols], use = 'complete.obs')
  
  
  ggplot(tab_sum) +
  geom_density(aes(x = WGT_VOL_Y), col ='yellow') +
  geom_density(aes(x = WGT_VOL_Q), col = 'blue') +
  geom_density(aes(x = VOL_Q), col = 'red') +
  geom_density(aes(x = VOL_Y), col = 'black')  +
  scale_x_continuous(trans = 'log')  
  
  
  write.csv(tab_sum,here::here('data_preparation_sampling',folder,"data_summary_by_siren.csv"))
  
  ############################# Plots of relevant distribution ###############################
  
  ################################### Plots of growth rates ##################################
  
  ################################# Yearly and quarterly GRs #################################
  
  gr_y_plot <- gather(gr_base[,1:26],YEAR,GR,`1994`:`2017`)
  gr_y_plot <- cbind(gr_y_plot,POOLTYPE = rep('Yearly logGR',nrow(gr_y_plot)))
  
  gr_q_plot <- gather(gr_base_q[,1:98],YEAR,GR,`1994.01`:`2017.04`) %>% 
               mutate(QUARTER = substr(YEAR,6,8)) %>% mutate(YEAR = substr(YEAR,1,4))
  gr_q_plot <- cbind(gr_q_plot,POOLTYPE = rep('Quarter-quarter logGR yearly-pooled',nrow(gr_q_plot)))
  
  gr_plot <- rbind(gr_y_plot, gr_q_plot[,c(1:4,6)])  
  
  rm(gr_y_plot)
  rm(gr_q_plot)
  
  fmt_decimals <- function(decimals=0)
  {
    function(x) formatC(x, format = 'e', digits = decimals)#signif(x, digits=decimals)
  }
  
  setDT(gr_plot)
  p <- gr_plot[YEAR %in% c(1994,2001,2010,2017)] %>%
       ggplot() +
       geom_density(aes(x=GR, group = YEAR, color = YEAR),n = 35) +
       scale_y_continuous(trans = 'log', limits = c(1e-7,1), labels = fmt_decimals(2)) +
       xlim(-10.8,10.8) +
       facet_wrap(~POOLTYPE, ncol = 2) 
      
  pp <- ggplot_build(p)
  
  groupsub <- function(x)
  {
    ifelse(x==1,1994,ifelse(x==2,2001,ifelse(x==3,2010,ifelse(x==4,2017,x))))
  }
  
  pp$data[[1]]$group <- factor(groupsub(pp$data[[1]]$group))
  
  colnames(pp$data[[1]])[colnames(pp$data[[1]]) == 'group'] <- 'YEAR'
  
  pacman::p_load(normalp)
  
  ss <-(10.8+10.8)/50
  
  data_base <- as.data.frame(cbind(Laplace = dnormp(seq(-10.8,10.8,ss),0,1,1),
                                   Gaussian = dnormp(seq(-10.8,10.8,ss),0,1,2), 
                                   x = seq(-10.8,10.8,ss)))
  
  data_base <- gather(data_base,type,y,1:2)
  
  
  pdf(here::here('data_preparation_sampling',folder,'logGR_density.pdf'))
  p +   
    geom_point(data = pp$data[[1]],aes(x = x, y = density, group = YEAR, col = YEAR, shape = YEAR), size = 1.5) +
    geom_line(data = data_base,aes(x = x, y = y, linetype = type, group = type), alpha = 0.5) +
    scale_color_hue(l=40, c=65) +
    theme_minimal() +
    scale_shape_manual(values=c(15, 16, 17, 18)) +
    guides(linetype=guide_legend(title="Benchmarks")) +
    labs(y = 'Binned density', x = 'Log-Growth Rates', title = '') 
  
  print(p)
  dev.off()
  
  
  ############################## Number of countries per year #################################
  
  plot_num <- num_sum_y[,c('SIREN','PYOD',paste0('NUM_',1993:2017))]
  
  colnames(plot_num)[3:ncol(plot_num)] <- 1993:2017
  
  plot_num <- plot_num[,-c('PYOD')]
  plot_num <- unique(plot_num)
  plot_num_long <- gather(plot_num,YEAR,NUMBER,`1993`:`2017`)
  
  
  abc <- plot_num_long %>% filter(YEAR %in% c(1993,2000,2009,2017)) 
  setDT(abc)  
  
  pdf(here::here('data_preparation_sampling',folder,'log_log_num_destinations.pdf'))
  p <- abc[,'N_FIRMS' := .N, keyby = c('NUMBER','YEAR')] %>% filter(NUMBER !=0) %>%
        select('YEAR','NUMBER','N_FIRMS')                 %>%
        unique()                                          %>%
        ggplot(aes(x = NUMBER, y = N_FIRMS))+
        geom_point(aes(shape = YEAR, col = YEAR), size = 2, alpha = 0.7) +
        labs(y = 'Num. of firms', x = 'Num. of destinations', title = '') +
        scale_color_hue(l=40, c=35) +
        theme_minimal() +
        scale_shape_manual(values=c(15, 16, 17, 18))+
        guides(col=guide_legend(title="Years")) +
        guides(shape=guide_legend(title="Years")) +  
        scale_y_continuous(trans = 'log10')  +
        scale_x_continuous(trans = 'log2')  +
        facet_wrap(~YEAR, ncol = 2) 
     
  print(p)
 
  dev.off()
  
  
  ####################################### Shannon index ########################################
  
  sh_sum_y <- sh_sum_y[,c('SIREN','PYOD',paste0('SH_',cols_y)), with = F][,paste0('S_IN',cols_y) := lapply(.SD,function(x)S_in_mod(x)),
                                                                          keyby = c('SIREN'),
                                                                          .SDcols = paste0('SH_',cols_y)]
  
  plot_shannon <- sh_sum_y[,c('SIREN','PYOD',paste0('S_IN',1993:2017))]
  colnames(plot_shannon)[3:ncol(plot_shannon)] <- 1993:2017
  
  plot_shannon <- plot_shannon[,-c('PYOD')]
  plot_shannon <- unique(plot_shannon)
  plot_shannon_long <- gather(plot_shannon,YEAR,NUMBER,`1993`:`2017`)
  
  
  pdf(here::here('data_preparation_sampling',folder,'shannon_index_density_4.pdf'))
 p <- plot_shannon_long %>% filter(YEAR %in% c(1993,2000,2009,2017)) %>% filter(NUMBER !=0) %>%
    ggplot(aes(x = NUMBER))+
    geom_density(aes(col = YEAR), alpha = 0, n = 100) +
    labs(y = 'Density', x = 'Shannon Index', title = '') +
    scale_color_hue(l=25, c=110) +
    theme_minimal() +
    guides(col=guide_legend(title="Years")) +
    guides(fill=guide_legend(title="Years")) +
    facet_wrap(~YEAR, ncol = 2)  

  print(p)
 
  dev.off()
  
  ####################################### Herfindal index ########################################
  
  sh_sum_y <- sh_sum_y[,c('SIREN','PYOD',paste0('SH_',cols_y)), with = F][,paste0('H_IN',cols_y) := lapply(.SD,function(x)H_in_mod(x)),
                                                                          keyby = c('SIREN'),
                                                                          .SDcols = paste0('SH_',cols_y)]
  
  
  plot_herfindal <- sh_sum_y[,c('SIREN','PYOD',paste0('H_IN',1993:2017))]
  colnames(plot_herfindal)[3:ncol(plot_herfindal)] <- 1993:2017
  
  plot_herfindal <- plot_herfindal[,-c('PYOD')]
  plot_herfindal <- unique(plot_herfindal)
  plot_herfindal_long <- gather(plot_herfindal,YEAR,NUMBER,`1993`:`2017`)
  
  pdf(here::here('data_preparation_sampling', folder,'herfindahl_index_density_4.pdf'))
p <-   plot_herfindal_long %>% filter(YEAR %in% c(1993,2000,2009,2017)) %>% filter(NUMBER != 1) %>%
    ggplot(aes(x = NUMBER))+
    geom_density(aes(fill = YEAR, col = YEAR), alpha = 0.8, n = 200) +
    labs(y = 'Density', x = 'Herfindahl Index', title = '') +
    scale_color_hue(l=40, c=25) +
    scale_fill_hue(l=40, c=25) +  
    theme_minimal() +
    guides(col=guide_legend(title="Years")) +
    guides(fill=guide_legend(title="Years")) +
    facet_wrap(~YEAR, ncol = 2)  
  
  print(p)

  dev.off()
  
  pdf(here::here('data_preparation_sampling', folder,'herfindahl_index_density_all.pdf'))
  p <-  plot_herfindal_long %>% filter(YEAR %in% c(1993,2000,2009,2017)) %>% filter(NUMBER != 1) %>%
    ggplot(aes(x = NUMBER))+
    geom_density(aes( col = YEAR), alpha = 0.8, n = 200) +
    labs(y = 'Density', x = 'Herfindahl Index', title = '') +
    scale_color_hue(l=25, c=110) +
    theme_minimal() +
    guides(col=guide_legend(title="Years")) 
  print(p)
  dev.off()
  
  ####################################### Herfindal index ########################################
  
  sh_sum_y <- sh_sum_y[,c('SIREN','PYOD',paste0('SH_',cols_y)), with = F][,paste0('DIV_IN',cols_y) := lapply(.SD,function(x)DIV_in_mod(x)),
                                                                          keyby = c('SIREN'),
                                                                          .SDcols = paste0('SH_',cols_y)]
  
  
  plot_herfindal <- sh_sum_y[,c('SIREN','PYOD',paste0('DIV_IN',1993:2017))]
  colnames(plot_herfindal)[3:ncol(plot_herfindal)] <- 1993:2017
  
  plot_herfindal <- plot_herfindal[,-c('PYOD')]
  plot_herfindal <- unique(plot_herfindal)
  plot_herfindal_long <- gather(plot_herfindal,YEAR,NUMBER,`1993`:`2017`)
  
  pdf(here::here('data_preparation_sampling', folder,'herfindahl_index_density_4.pdf'))
  p <-   plot_herfindal_long %>% filter(YEAR %in% c(1993,2000,2009,2017)) %>% filter(NUMBER != 1) %>%
    ggplot(aes(x = NUMBER))+
    geom_density(aes(fill = YEAR, col = YEAR), alpha = 0.8, n = 200) +
    labs(y = 'Density', x = 'Diversification Index', title = '') +
    scale_color_hue(l=40, c=25) +
    scale_fill_hue(l=40, c=25) +  
    theme_minimal() +
    guides(col=guide_legend(title="Years")) +
    guides(fill=guide_legend(title="Years")) +
    facet_wrap(~YEAR, ncol = 2)  
  
  print(p)
  
  dev.off()
  
  pdf(here::here('data_preparation_sampling', folder,'herfindahl_index_density_all.pdf'))
  p <-  plot_herfindal_long %>% filter(YEAR %in% c(1993,2000,2009,2017)) %>% filter(NUMBER != 1) %>%
    ggplot(aes(x = NUMBER))+
    geom_density(aes( col = YEAR), alpha = 0.8, n = 200) +
    labs(y = 'Density', x = 'Diversification Index', title = '') +
    scale_color_hue(l=25, c=110) +
    theme_minimal() +
    guides(col=guide_legend(title="Years")) 
  print(p)
  dev.off()
  
  ####################################### Share of top export ##################################
  
  sh_sum_y <- sh_sum_y[,c('SIREN','PYOD',paste0('SH_',cols_y)), with = F][,paste0('S_TOP',cols_y) := lapply(.SD,function(x)max(x)),
                                                                          keyby = c('SIREN'),
                                                                          .SDcols = paste0('SH_',cols_y)]
  
  
  plot_s_top <- sh_sum_y[,c('SIREN','PYOD',paste0('S_TOP',1993:2017))]
  colnames(plot_s_top)[3:ncol(plot_s_top)] <- 1993:2017
  
  plot_s_top <- plot_s_top[,-c('PYOD')]
  plot_s_top <- unique(plot_s_top)
  plot_s_top_long <- gather(plot_s_top,YEAR,NUMBER,`1993`:`2017`)
  
  pdf(here::here('data_preparation_sampling',folder,'s_top_index_density_4.pdf'))
  p <-  plot_s_top_long %>% filter(YEAR %in% c(1993,2000,2009,2017)) %>% filter(NUMBER != 0) %>%
        ggplot(aes(x = NUMBER))+
        geom_density(aes(fill = YEAR, col = YEAR), alpha = 0.8, n = 200) +
        labs(y = 'Density', x = 'Share of the most important destination', title = '') +
        scale_color_hue(l=40, c=25) +
        scale_fill_hue(l=40, c=25) +  
        theme_minimal() +
        guides(col=guide_legend(title="Years")) +
        guides(fill=guide_legend(title="Years")) +
        facet_wrap(~YEAR, ncol = 2)  

  print(p)

  dev.off()
  
  
  pdf(here::here('data_preparation_sampling',folder,'s_top_index_density_all.pdf'))
  p <-  plot_s_top_long %>% filter(YEAR %in% c(1993,2000,2009,2017)) %>% filter(NUMBER != 0) %>%
        ggplot(aes(x = NUMBER))+
        geom_density(aes( col = YEAR), alpha = 0.8, n = 200) +
        labs(y = 'Density', x = 'Share of the most important destination', title = '') +
        scale_color_hue(l=25, c=110) +
        theme_minimal() +
        guides(col=guide_legend(title="Years")) 
     
  print(p)
   
  dev.off()
}

statistics_and_volatility <- function(tab_sum, folder = 'data_summary')
  {
  
  
  ##### Create the folders where to store the outcomes (if not existing) #####
  dir.create(here::here('data_preparation_sampling'), showWarnings = FALSE)
  dir.create(here::here('data_preparation_sampling', folder), showWarnings = FALSE)
  
  ##### Auxiliary packages ####     
  
  library(ggridges)
  library(psych)
  library(xtable)
  
  
  ######################### Descriptive statistics on different time units ##################
  
  tab_cols <- colnames(tab_sum)[2:ncol(tab_sum)]
  tab_des <- tab_sum[,lapply(.SD,as.numeric) ,.SDcols = tab_cols]

  whole_ts <- tab_des[,c('WGT_VOL_Y','WGT_VOL_Q','NUM','LOG_NUM','S_TOP','H_IN','DIV_IN','S_IN')]
  whole_ts <- na.omit(whole_ts)

  
  yearly_basis <- tab_des[,c('WGT_VOL_Y','WGT_VOL_Q','MEAN_NUM_Y','MEDIAN_NUM_Y','S_TOP_MEAN_Y','H_IN_MEAN_Y','DIV_IN_MEAN_Y','S_IN_MEAN_Y')]
  yearly_basis <- na.omit(yearly_basis)
  
  quarterly_basis <- tab_des[,c('WGT_VOL_Y','WGT_VOL_Q','MEAN_NUM_Q','MEDIAN_NUM_Q','S_TOP_MEAN_Q','H_IN_MEAN_Q','DIV_IN_MEAN_Q','S_IN_MEAN_Q')]
  quarterly_basis <- na.omit(quarterly_basis)
  

  whole_sum <- psych::describe(whole_ts)[,c(2:5,8:9,11:12)]
  colnames(whole_sum) <- c('N. obs.','Mean','SD','Median','Min','Max','Skewness','Kurtosis')
  rownames(whole_sum) <- c("$\\sigma_y(g)$",'$\\sigma_q(g)$','# of dest.','$\\log(#)$','Top share','Herf. index','Div. index','Shannon index')
  
  
  yearly_sum <- psych::describe(yearly_basis)[,c(2:5,8:9,11:12)]
  colnames(yearly_sum) <- c('N. obs.','Mean','SD','Median','Min','Max','Skewness','Kurtosis')
  rownames(yearly_sum) <- c("$\\sigma_y(g)$",'$\\sigma_q(g)$','Mean # of dest.','$Median # of dest.$','Mean top sh.','Mean Herf. index','Mean div. index','Mean Shan. index')
  
  
  quarterly_sum <- psych::describe(quarterly_basis)[,c(2:5,8:9,11:12)]
  colnames(quarterly_sum) <- c('N. obs.','Mean','SD','Median','Min','Max','Skewness','Kurtosis')
  rownames(quarterly_sum) <- c('$\\sigma_y(g)$','$\\sigma_q(g)$','Mean # of dest.','$Median # of dest.$','Mean top sh.','Mean Herf. index','Mean div. index','Mean Shan. index')
  
  
  print(xtable(whole_sum),file=here::here('data_preparation_sampling',folder,"time_span_summary.txt"))
  print(xtable(yearly_sum),file=here::here('data_preparation_sampling',folder,"yearly_summary.txt"))
  print(xtable(quarterly_sum),file=here::here('data_preparation_sampling',folder,"quarterly_summary.txt"))
  
  
  upper_clean <- function(x)
  {
    for(i in 1:(nrow(x)-1))
    {   
      for(j in (i+1):ncol(x)){x[i,j] <- ''}  
    }
    return(x)
  }
  
  cor_mat_whole <- upper_clean(format(round(cor(whole_ts, use = 'complete.obs'),4), nsmall = 4))
  colnames(cor_mat_whole) <- c("\\sigma_y(g)",'\\sigma_q(g)','# of dest.','$\\log(#)$','Top share','Herf. index','Div. index','Shannon index')
  rownames(cor_mat_whole) <- c("\\sigma_y(g)",'\\sigma_q(g)','# of dest.','$\\log(#)$','Top share','Herf. index','Div. index','Shannon index')
  
  cor_mat_y <- upper_clean(format(round(cor(yearly_basis, use = 'complete.obs'),4), nsmall = 4))
  colnames(cor_mat_y) <- c("\\sigma_y(g)",'\\sigma_q(g)','Mean # of dest.','Median # of dest.','Mean top share','Mean Herf. index','Mean div. index','Mean Shannon index')
  rownames(cor_mat_y) <- c("\\sigma_y(g)",'\\sigma_q(g)','Mean # of dest.','Median # of dest.','Mean top share','Mean Herf. index','Mean div. index','Mean Shannon index')
  
  cor_mat_q <- upper_clean(format(round(cor(quarterly_basis, use = 'complete.obs'),4), nsmall = 4))
  colnames(cor_mat_q) <- c("\\sigma_y(g)",'\\sigma_q(g)','Mean # of dest.','Median # of dest.','Mean top share','Mean Herf. index','Mean div. index','Mean Shannon index')
  rownames(cor_mat_q) <- c("\\sigma_y(g)",'\\sigma_q(g)','Mean # of dest.','Median # of dest.','Mean top share','Mean Herf. index','Mena div. index','Mean Shannon index')
  
  print(xtable(cor_mat_whole),file=here::here('data_preparation_sampling',folder,"time_span_correlations.txt"))
  print(xtable(cor_mat_y),file=here::here('data_preparation_sampling',folder,"yearly_correlations.txt"))
  print(xtable(cor_mat_q),file=here::here('data_preparation_sampling',folder,"quarterly_correlations.txt"))
  
  
  #### PLotting auxiliary functions ####   
  
  fmt_decimals2 <- function(decimals=0)
  {
    function(x) signif(x, digits=decimals)
  }
  binning <- function(x)
  {
    
    x <- as.numeric(x)
    qv <- quantile(x,(1:5/5), na.rm = T)
    x[x>qv[4]] <- '5th quintile'
    x[x<=qv[1]] <- '1st quintile'
    x[x>qv[1] & x<=qv[2]] <- '2nd quintile'
    x[x>qv[2] & x<=qv[3]] <- '3rd quintile'
    x[x>qv[3] & x<=qv[4]] <- '4th quintile'
    
    return(factor(x, levels = c('1st quintile','2nd quintile', '3rd quintile',
                                '4th quintile', '5th quintile')))
  }
  binning_auto <- function(x)
  {
    
    x <- as.numeric(x)
    qv <- quantile(x,(1:9/10), na.rm = T)
    qv <- unique(qv)
    
    aux_df <- data.frame(x)
    setDT(aux_df)
    aux_df[,q := cut(x,breaks = c(-Inf,qv, Inf), na.rm = T, include.lowest = T)] 
    
    
    return(aux_df$q)
  }
  
  ################### Volatility versus concentration and size measures ######################

  data_sum <- tab_sum[,c('NUM','LOG_NUM','WGT_VOL_Y','WGT_VOL_Q')] %>%
              gather(vol_type,vol_val,3:4) %>%
              mutate(BIN_LOG_NUM = binning_auto(LOG_NUM))
  setDT(data_sum)      
  
  #### PLotting distributions logNUM vs vol
  pdf(here::here('data_preparation_sampling',folder,'logN_vs_vol_density.pdf'))
  data_sum$BIN_LOG_NUM
  p <- data_sum[!is.na(BIN_LOG_NUM)]  %>%  
      ggplot(aes(x = 0.01 + sqrt(vol_val), y = BIN_LOG_NUM)) +
      scale_x_continuous(labels = fmt_decimals2(4)) + 
      stat_density_ridges(aes(fill = vol_type, color = vol_type), alpha = 0.2, scale = 1.2,
                          quantiles = 10, quantile_lines = TRUE) +
      geom_boxplot(aes(fill = vol_type), alpha = 0.5, width = 0.06, outlier.shape = NA) + 
      scale_fill_hue(l=40, c=35,
                     labels = c("Volatility on quarterly basis","Volatility on yearly basis"))  +
      scale_color_hue(l=35, c=35, 
                      labels = c("Volatility on quarterly basis","Volatility on yearly basis")) + 
      labs(y = '', x = '', title = '') +  
      guides(fill=guide_legend(title=""),
             col =guide_legend(title="") ) +
      theme_minimal()  +
      theme(legend.direction = 'horizontal',
            legend.position = 'top')  
  print(p)
  dev.off()
  
  
  #### PLotting distributions herfindal vs vol
  

  data_sum_h <- tab_sum[,c('H_IN','WGT_VOL_Y','WGT_VOL_Q')] %>%
                gather(vol_type,vol_val,2:3) %>%
                mutate(BIN_H = binning_auto(H_IN))
  setDT(data_sum_h)

  #### PLotting distributions Hin vs vol
  
  pdf(here::here('data_preparation_sampling',folder,'Hin_vs_vol_density.pdf'))
  p <-  data_sum_h[!is.na(BIN_H)]  %>%  
        ggplot(aes(x = 0.01 + sqrt(vol_val), y = BIN_H)) +
        scale_x_continuous(labels = fmt_decimals2(4)) + 
        stat_density_ridges(aes(fill = vol_type, color = vol_type), alpha = 0.2, scale = 1,
                            quantiles = 4, quantile_lines = TRUE) +
        geom_boxplot(aes(fill = vol_type), alpha = 0.5, width = 0.06, outlier.shape = NA) + 
        scale_fill_hue(l=40, c=35,
                       labels = c("Volatility on quarterly basis","Volatility on yearly basis"))  +
        scale_color_hue(l=35, c=35, 
                        labels = c("Volatility on quarterly basis","Volatility on yearly basis")) + 
        labs(y = '', x = '', title = '') +  
        guides(fill=guide_legend(title=""),
               col =guide_legend(title="") ) +
        theme_minimal()  +
        theme(legend.direction = 'horizontal',
              legend.position = 'top')  
     
  print(p)
  dev.off()
  
  
  #### PLotting distributions diversification vs vol
  
  
  data_sum_h <- tab_sum[,c('DIV_IN','WGT_VOL_Y','WGT_VOL_Q')] %>%
                gather(vol_type,vol_val,2:3) %>%
                mutate(BIN_DIV = binning_auto(DIV_IN))
              setDT(data_sum_h)
  
  #### plot
  
  pdf(here::here('data_preparation_sampling',folder,'Div_index_vs_vol_density.pdf'))
  p <-  data_sum_h[!is.na(BIN_DIV)]  %>%  
    ggplot(aes(x = 0.01 + sqrt(vol_val), y = BIN_DIV)) +
    scale_x_continuous(labels = fmt_decimals2(4)) + 
    stat_density_ridges(aes(fill = vol_type, color = vol_type), alpha = 0.2, scale = 1,
                        quantiles = 4, quantile_lines = TRUE) +
    geom_boxplot(aes(fill = vol_type), alpha = 0.5, width = 0.06, outlier.shape = NA) + 
    scale_fill_hue(l=40, c=35,
                   labels = c("Volatility on quarterly basis","Volatility on yearly basis"))  +
    scale_color_hue(l=35, c=35, 
                    labels = c("Volatility on quarterly basis","Volatility on yearly basis")) + 
    labs(y = '', x = '', title = '') +  
    guides(fill=guide_legend(title=""),
           col =guide_legend(title="") ) +
    theme_minimal()  +
    theme(legend.direction = 'horizontal',
          legend.position = 'top')  
  
  print(p)
  dev.off()
  
  #### PLotting distributions diversification index vs vol (yearly)
  
  data_sum_h <- tab_sum[,c('DIV_IN_MEAN_Y','WGT_VOL_Y','WGT_VOL_Q')] %>%
                gather(vol_type,vol_val,2:3) %>%
                mutate(BIN_DIV_Y = binning_auto(DIV_IN_MEAN_Y))
  
  setDT(data_sum_h)
  
  
  
  pdf(here::here('data_preparation_sampling',folder,'Div_index_yearly_mean_vs_vol_density.pdf'))
  p <-  data_sum_h[!is.na(BIN_DIV_Y)]  %>%  
    ggplot(aes(x = 0.01 + sqrt(vol_val), y = BIN_DIV_Y)) +
    scale_x_continuous(labels = fmt_decimals2(4)) + 
    stat_density_ridges(aes(fill = vol_type, color = vol_type), alpha = 0.2, scale = 1,
                        quantiles = 4, quantile_lines = TRUE) +
    geom_boxplot(aes(fill = vol_type), alpha = 0.5, width = 0.06, outlier.shape = NA) + 
    scale_fill_hue(l=40, c=35,
                   labels = c("Volatility on quarterly basis","Volatility on yearly basis"))  +
    scale_color_hue(l=35, c=35, 
                    labels = c("Volatility on quarterly basis","Volatility on yearly basis")) + 
    labs(y = '', x = '', title = '') +  
    guides(fill=guide_legend(title=""),
           col =guide_legend(title="") ) +
    theme_minimal()  +
    theme(legend.direction = 'horizontal',
          legend.position = 'top')  
  
  print(p)
  dev.off()
  
  #### PLotting distributions diversification index vs vol (quarterly)
  
  data_sum_h <- tab_sum[,c('DIV_IN_MEAN_Q','WGT_VOL_Y','WGT_VOL_Q')] %>%
    gather(vol_type,vol_val,2:3) %>%
    mutate(BIN_DIV_Q = binning_auto(DIV_IN_MEAN_Q))
  
  setDT(data_sum_h)
  
  
  
  pdf(here::here('data_preparation_sampling',folder,'Div_index_quarterly_mean_vs_vol_density.pdf'))
  p <-  data_sum_h[!is.na(BIN_DIV_Q)]  %>%  
    ggplot(aes(x = 0.01 + sqrt(vol_val), y = BIN_DIV_Q)) +
    scale_x_continuous(labels = fmt_decimals2(4)) + 
    stat_density_ridges(aes(fill = vol_type, color = vol_type), alpha = 0.2, scale = 1,
                        quantiles = 4, quantile_lines = TRUE) +
    geom_boxplot(aes(fill = vol_type), alpha = 0.5, width = 0.06, outlier.shape = NA) + 
    scale_fill_hue(l=40, c=35,
                   labels = c("Volatility on quarterly basis","Volatility on yearly basis"))  +
    scale_color_hue(l=35, c=35, 
                    labels = c("Volatility on quarterly basis","Volatility on yearly basis")) + 
    labs(y = '', x = '', title = '') +  
    guides(fill=guide_legend(title=""),
           col =guide_legend(title="") ) +
    theme_minimal()  +
    theme(legend.direction = 'horizontal',
          legend.position = 'top')  
  
  print(p)
  dev.off()
  
  
  #################  Shannon index ##############################
  
  
  data_sum_s <- tab_sum[,c('S_IN','WGT_VOL_Y','WGT_VOL_Q')] %>%
                gather(vol_type,vol_val,2:3) %>% mutate(BIN_S = binning_auto(S_IN))
  
  setDT(data_sum_s)
  
  #### PLotting distributions SH index vs vol
  
  
  pdf(here::here('data_preparation_sampling', folder,'sh_vs_vol_density.pdf'))
  p <-  data_sum_s[!is.na(BIN_S)]  %>%  
        ggplot(aes(x = 0.01 + sqrt(vol_val), y = BIN_S)) +
        scale_x_continuous(labels = fmt_decimals2(4)) + 
        stat_density_ridges(aes(fill = vol_type, color = vol_type), alpha = 0.2, scale = 1,
                            quantiles = 4, quantile_lines = TRUE) +
        geom_boxplot(aes(fill = vol_type), alpha = 0.5, width = 0.06, outlier.shape = NA) + 
        scale_fill_hue(l=40, c=35,
                       labels = c("Volatility on quarterly basis","Volatility on yearly basis"))  +
        scale_color_hue(l=35, c=35, 
                        labels = c("Volatility on quarterly basis","Volatility on yearly basis")) + 
        labs(y = '', x = '', title = '') +  
        guides(fill=guide_legend(title=""),
               col =guide_legend(title="") ) +
        theme_minimal()  +
        theme(legend.direction = 'horizontal',
              legend.position = 'top')
     
   print(p)
   
   dev.off()
  
  
  #################  S_top index ##############################
  
  

  data_sum_stop <- tab_sum[,c('S_TOP','WGT_VOL_Y','WGT_VOL_Q')] %>%
                    gather(vol_type,vol_val,2:3) %>% mutate(BIN_STOP = binning_auto(S_TOP))
  
  setDT(data_sum_stop)
  
  #### PLotting distributions S_top vs vol
  
  pdf(here::here('data_preparation_sampling',folder,'s_top_vs_vol_density.pdf'))
  p <- data_sum_stop[!is.na(BIN_STOP)]  %>%  
      ggplot(aes(x = 0.01 + sqrt(vol_val), y = BIN_STOP)) +
      scale_x_continuous(labels = fmt_decimals2(4)) + 
      stat_density_ridges(aes(fill = vol_type, color = vol_type), alpha = 0.2, scale = 0.9,
                          quantiles = 4, quantile_lines = TRUE) +
      geom_boxplot(aes(fill = vol_type), alpha = 0.5, width = 0.06, outlier.shape = NA) + 
      scale_fill_hue(l=40, c=35,
                     labels = c("Volatility on quarterly basis","Volatility on yearly basis"))  +
      scale_color_hue(l=35, c=35, 
                      labels = c("Volatility on quarterly basis","Volatility on yearly basis")) + 
      labs(y = '', x = '', title = '') +  
      guides(fill=guide_legend(title=""),
             col =guide_legend(title="") ) +
      theme_minimal()  +
      theme(legend.direction = 'horizontal',
            legend.position = 'top')  
    print(p)
    dev.off()
 
}




