################################################################################
############ Main script for results analysis after model estimation ###########
################# Results for the firm volatility - Micro level ################
################################################################################

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

#### Auxiliary packages ####
pacman::p_load(here)
pacman::p_load(data.table)
pacman::p_load(tidyverse)
pacman::p_load(ggplot2)
pacman::p_load(ggridges)
pacman::p_load(gridExtra)
pacman::p_load(matrixStats)
pacman::p_load(stargazer)
pacman::p_load(xtable)



#### Auxiliary functions ####

f_one = function(DT) 
{
  for (i in names(DT))
    DT[!is.na(get(i)), (i):=1]
}
f_zero = function(DT) 
{
  for (i in names(DT))
    DT[is.na(get(i)), (i):=0]
}
f_nan = function(DT) 
{
  for (i in names(DT))
    DT[get(i)==0, (i):=NA]
}
fmt_decimals <- function(decimals=0) ## Auxiliary function to format the numbers in the plot axis
{
  function(x) formatC(x, format = 'e', digits = decimals)#signif(x, digits=decimals)
}    
fmt_decimals2 <- function(decimals=0)
{
  function(x) signif(x, digits=decimals)
}


bootstrap_CI <- function(x,q)
{
  avg <- mean(x, na.rm = T)
  vec_delta <- colMeans(sapply(1:100,function(y) sample(x,length(x), replace = T)))-avg 
  CI <- sort(vec_delta)[floor(q*100)]
  return(CI)     
}
ci_fun <- function(vec)
{
  alpha <- 0.01
  se <- sd(vec)/sqrt(length(vec))
  t <- qt((1-alpha)/2 + .5, length(vec)-1)
  CI = t*se
  return(CI)
}

## aggregation method as in DiGiovanni et al.
## variance calculations and related standard errors

const_weights_aggregate <- function(m_weights, m_gr)
{
  output_list <- list() ## define the list to store the output
  
  m_weights[is.na(m_weights)] <- 0
  m_gr <- m_gr - rowMeans(m_gr, na.rm = T)
  m_gr[is.na(m_gr)] <- 0
  counts <- m_gr
  counts[m_gr!= 0] <- 1
  
  m_weights <- m_weights[rowSums(counts) > 2,]
  m_gr <- m_gr[rowSums(counts) > 2,]
  counts <- counts[rowSums(counts) > 2,]
  
  counts_mat <- matrix(rep(rowSums(counts), ncol(counts)), ncol = ncol(counts))
  m_gr <- as.matrix(m_gr/sqrt(counts_mat - 1))
  
  var_mat <- t(m_weights)%*%m_gr
  
  output_list[['var_mat']] <- var_mat
  output_list[['var']] <- diag(var_mat%*%t(var_mat))
  output_list[['sd']] <- sqrt(output_list$var)
  
  TT <- dim(var_mat)[1] ## number of time steps
  q <- round(.75*(TT)^(1/3) - 1) ## approximate the
  
  dif_c <- (sqrt(TT-1)*(output_list$var)^2)-(1/sqrt(TT-1)*matrix(rep(output_list$var,TT), ncol = TT))
  S2 <- rowSums(dif_c^2)
  
  ## autorcorreltions
  
  S12 <- matrix(0,TT,1)
  for (k in 1:q)
  {  
    c <- 2*(1 - k/(q+1))
    dif_Tk <- (sqrt(TT-k-1)*(var_mat^2))-(1/sqrt(TT-k-1)*matrix(rep(output_list$var,TT), ncol = TT))
    dif_Tmk <- dif_Tk[,1:(ncol(dif_Tk)-k)]
    dif_Tpk <- dif_Tk[,(1+k):ncol(dif_Tk)]
    s12 <- rowSums(dif_Tmk*dif_Tpk)
    SS12 <- c*s12
    S12 <- S12 + SS12
  }
  
  
  se_var <- sqrt(abs(S2 + S12))/sqrt(TT)
  se_sd <- se_var/(2*sqrt(output_list$var))
  
  output_list[['se_var']] <- se_var
  output_list[['se_sd']] <- se_sd
  
  return(output_list)
} 


### Folder containing the results of the estimation procedure
folder_res <- 'results_samples_26_5k' ### folder_res <- 'results_samples_68_100k'

## Check if the folder with the results is correctly specified
folds_list <- list.files(here::here('results_analysis','Inputs'))
if(!sum(grepl(folder_res,folds_list)))
{
  cat('\nThere is no folder named: ', folder_res)
  folder_res <- folds_list[1]  
  cat('\nUsing the first available folder instead: ',folder_res)
}


### Subfolder containing the samples
folder_samp <- 'samples'

#### Cycle setup  ####

siren_vol_list <- list()

H <- 0 #### number of dataset subsamples (0 means that we are workking with the full dataset only)

#### Cycling over the estimates from each sample of the dataset (including also the complete dataset i = 0)  ####

### Output storage lists
fac_list <- list()

agg_cor_list <- list()
agg_vol_list <- list()

dest_vol_list <- list()

vol_comparison_ts_list <- list()
vol_sample_ts_list <- list()
matrix_list <- list()

for(i in 0)
{
  cat('\nStarting the step',i,'of the cycle',ifelse(i==0,'(at step 0 the complete dataset is used)\n','\n'))
  
  ## Import firm weights
  cat('   Importing the firm\'s size-based weights...\n')
  w_firm <- read_csv(here::here('results_analysis','Inputs',folder_res,folder_samp,paste0('data_sample_w_firm',i,'.csv')))
  setDT(w_firm)
  w_firm <- w_firm[,-1]
  
  
  ############################# Firms' export features #################################
  
  cat('   Calculating diversification measures at the firm-level...\n')
  
  ##### Diversification measures
   H_in <- function(x) ## given a vector of shares calculate the Herfindal index
   {
     ifelse(sum(!is.na(x))==0,-1,sum(x^2, na.rm = T))
   }
   
   S_top <- function(x)
   {
     ifelse(sum(!is.na(x))==0,-1,max(x, na.rm = T))
   }
   
   Var_in <- function(x) ### auxiliary function evaluating coherence of the number of destinations
   {
     aa <- which(colSums(x)>0)
     n_min <- min(aa)
     n_max <- max(aa)
     
     na_span <- sum(!colSums(x[,n_min:n_max])>0)
     span <- n_max-n_min+1
     
     x <- as.matrix(x[,n_min:n_max])
     
     span_N <- mean(colSums(x))
     N <- nrow(x)
     
     y <- x%*%t(x)
     ind <- (sum(y)+sum(diag(y)))/2          
     
     return(c(N,span_N,as.integer(span),as.integer(na_span),ind/((n_max-n_min)*dim(x)[1])))
   }

   ts_cols <- colnames(w_firm)[3:ncol(w_firm)]
   cols <- colnames(w_firm)[1:ncol(w_firm)]
  # 
  # 
  # ####### Herf index in various flavours
   cat('     Calculating the HHI at the firm-level...\n')
   h_in <- w_firm[, lapply(.SD, H_in), keyby = SIREN, .SDcols = ts_cols] ### calculate herf ind column by column    h_in[h_in==-1] <- NA
   h_in[h_in==-1] <- NA
   h_in <- cbind(h_in[,1],H0 = h_in[,2], HI = rowMeans(h_in[,2:5], na.rm = T), HQ = rowMeans(h_in[,2:ncol(h_in)], na.rm =T),
                 C0 = 1/h_in[,2], CI = rowMeans(1/h_in[,2:5], na.rm = T), CQ = rowMeans(1/h_in[,2:ncol(h_in)], na.rm = T))
   colnames(h_in)[c(2,5)] <- c('H0','C0')
   
   cat('     Calculating the top share at the firm-level...\n')
   s_top <- w_firm[, lapply(.SD, S_top), keyby = SIREN, .SDcols = ts_cols] ### calculate top share column by column
   s_top[s_top==-1] <- NA
   s_top <- cbind(s_top[,1],ST0 =s_top[,2], STI = rowMeans(s_top[,2:5], na.rm = T), STQ = rowMeans(s_top[,2:ncol(h_in)], na.rm =T))
   colnames(s_top)[2] <- c('ST0')
   
   
  # cat('     Calculating diversification measures based on the number of destinations at the firm-level...\n')
  # 
  # x  <- w_firm  
  # x[is.na(x)] <- 0  
  # for(s in unique(w_firm$SIREN))
  # {
  #   y <- x[SIREN == s,-c(1,2)]
  #   y[y>0] <- 1
  #   
  #   abc <- Var_in(y)
  #   h_in[SIREN == s,c('NUM')] <- abc[1]
  #   h_in[SIREN == s,c('NUMQ')] <- abc[2]
  #   h_in[SIREN == s,c('TLIFE')] <- abc[3]
  #   h_in[SIREN == s,c('NALIFE')] <- abc[4]
  #   h_in[SIREN == s,c('DDIQ')] <- abc[5]
  # }  
  # 
   setkey(s_top,SIREN)
   setkey(h_in, SIREN)
   
   h_in <- s_top[h_in]
  
  
  ######################### Imports and components definition and calculations #########################
  
  ### Import factors
  cat('   Importing the factors\'s time series...\n')
  
  fac <- read_csv(here::here('results_analysis','Inputs',folder_res,paste0('results_',i,'_F.csv')),
                  show_col_types = FALSE)
  fac_tab <- gather(cbind(TIME = 1:nrow(fac),fac),FACTOR,VALUE,2:(ncol(fac)+1))
  fac_tab <- cbind(SIM = rep(i,nrow(fac_tab)), fac_tab)
  
  
  ## Import factor loadings
  cat('   Importing the loadings...\n')
  
  load <- setDT(read_csv(here::here('results_analysis','Inputs',folder_res,paste0('results_',i,'_C.csv')),
                         show_col_types = FALSE))
  
  ## Calculating the common component
  com_comp <- as.matrix(load[,3:ncol(load)])%*%t(as.matrix(fac)) ## common and dest components in one
  
  ## Import means and variances
  cat('   Importing means and variances of the time series...\n')
  
  sig <- read_csv(here::here('results_analysis','Inputs',folder_res,paste0('results_',i,'_W.csv')),
                  show_col_types = FALSE) ### vector of standard deviations
  means <- read_csv(here::here('results_analysis','Inputs',folder_res,paste0('results_',i,'_M.csv')),
                    show_col_types = FALSE) ### vector of means
  
  
  v_sig <- as.matrix(sig[,3])
  v_means <-  as.matrix(means[,3])
  

  weights_matrix <- as.matrix(w_firm[,3:(ncol(w_firm)-4)])
  
  ## Matrix with NA pattern
  w_na <- as.matrix(w_firm[,3:(ncol(w_firm)-4)])
  w_na[!is.na(w_na)] <- 1
  
  ## Matrix with time-series specific standard deviations and means
  ## (the NA pattern of the original matrix is imposed already at this stage)
  sd_mat_na <- matrix(rep(v_sig,ncol(weights_matrix)), ncol = ncol(weights_matrix))*w_na  
  mean_mat_na <- matrix(rep(v_means,ncol(weights_matrix)), ncol = ncol(weights_matrix))*w_na 
  
  
  ## Import flow growth rates
  cat('   Importing the growth rates dataset...\n')
  gr_mat <- setDT(read_csv(here::here('results_analysis','Inputs',folder_res,folder_samp,paste0('data_sample_gr',i,'.csv')),
                           show_col_types = FALSE,
                           name_repair = "minimal"))
  gr_mat <- gr_mat[,-1]
  
  
  cat('   Computing the weigthed decomposition for each flow...\n')
  ### Demeaning the actual growth rates (after demeaning)
  gr_actual_mat <- gr_mat[,3:ncol(gr_mat)] - rowMeans(gr_mat[,3:ncol(gr_mat)], na.rm = T)
  
  ## common (global+dest) component of the growth rates rescaled including the variance of the original series
  common_mat_na <-  sd_mat_na*as.matrix(load[,3:ncol(load)])%*%t(as.matrix(fac[,1:ncol(fac)]))
  
  ## common global component of the growth rates rescaled including the variance of the original series
  common_global_mat_na <- sd_mat_na*as.matrix(load[,3])%*%t(as.matrix(fac[,1])) 

  ## common destination component of the growth rates rescaled including the variance of the original series
  common_destination_mat_na <- sd_mat_na*as.matrix(load[,4:ncol(load)])%*%t(as.matrix(fac[,2:ncol(fac)])) 

  ## idiosyncratic component of the growth rates calculated by difference from the actual growth rates and rescaled common component
  idio_mat <- as.matrix(gr_actual_mat - common_mat_na) 
  
  
  cat('   Aggregating at the firm level...\n')
  
  ## Actual weighted growth rates matrix
  siren_wgt <- cbind(gr_mat[,1:2],gr_actual_mat*weights_matrix)
  agg_cols <- colnames(siren_wgt)[-c(1,2)]
  
  siren_series <- siren_wgt[, lapply(.SD, function(x) sum(x, na.rm = T)),
                            , keyby = c('SIREN'), .SDcols = agg_cols]
  
  siren_vol <- cbind(siren_series[,1],VOL = rowVars(as.matrix(siren_series[,-1]), na.rm = T))
  
  ## Component matrices
  siren_wgt_glob <- cbind(gr_mat[,1:2],common_global_mat_na*weights_matrix)
  siren_wgt_cnt <- cbind(gr_mat[,1:2],common_destination_mat_na*weights_matrix)
  siren_wgt_idio <- cbind(gr_mat[,1:2],idio_mat*weights_matrix)
  colnames(siren_wgt_glob) <- colnames(siren_wgt_cnt) <- colnames(siren_wgt_idio) <- c('SIREN','PYOD',agg_cols)
  
  ## Component series
  siren_glob_series <- siren_wgt_glob[, lapply(.SD, function(x) sum(x, na.rm = T)),
                                    , keyby = c('SIREN'), .SDcols = agg_cols]
  siren_cnt_series <- siren_wgt_cnt[, lapply(.SD, function(x) sum(x, na.rm = T)),
                                    , keyby = c('SIREN'), .SDcols = agg_cols]
  siren_idio_series <- siren_wgt_idio[, lapply(.SD, function(x) sum(x, na.rm = T)),
                                      , keyby = c('SIREN'), .SDcols = agg_cols]
  ## Muted series
  siren_mut_idio <- cbind(siren_glob_series[,1],siren_glob_series[,-1] + siren_cnt_series[,-1])
  siren_mut_glob <- cbind(siren_glob_series[,1],siren_idio_series[,-1] + siren_cnt_series[,-1])
  siren_mut_cnt <- cbind(siren_glob_series[,1],siren_glob_series[,-1] + siren_idio_series[,-1])
  
  
  cat('   Computing the volatilities at the firm level...\n')
  ## Component volatilities
  siren_glob_vol <- cbind(siren_glob_series[,1], GLOB_VOL = rowVars(as.matrix(siren_glob_series[,-1]), na.rm = T))
  siren_cnt_vol <- cbind(siren_cnt_series[,1], CNT_VOL = rowVars(as.matrix(siren_cnt_series[,-1]), na.rm = T))
  siren_idio_vol <- cbind(siren_idio_series[,1], IDIO_VOL = rowVars(as.matrix(siren_idio_series[,-1]), na.rm = T))
  
  ## Muted volatilities
  siren_mut_idio_vol <- cbind(siren_mut_idio[,1], MUT_IDIO_VOL = rowVars(as.matrix(siren_mut_idio[,-1]), na.rm = T))
  siren_mut_glob_vol <- cbind(siren_mut_glob[,1], MUT_GLOB_VOL = rowVars(as.matrix(siren_mut_glob[,-1]), na.rm = T))
  siren_mut_cnt_vol <- cbind(siren_mut_cnt[,1], MUT_CNT_VOL = rowVars(as.matrix(siren_mut_cnt[,-1]), na.rm = T))
  
  
  cat('   Building the firm-level dataset...\n')
  
  setkey(siren_vol, SIREN)
  setkey(siren_glob_vol, SIREN)
  setkey(siren_cnt_vol, SIREN)
  setkey(siren_idio_vol, SIREN)
  setkey(siren_mut_glob_vol, SIREN)
  setkey(siren_mut_cnt_vol, SIREN)
  setkey(siren_mut_idio_vol, SIREN)
  
  siren_vol <- merge(siren_vol,siren_glob_vol, all = T)
  siren_vol <- merge(siren_vol,siren_cnt_vol, all = T)
  siren_vol <- merge(siren_vol,siren_idio_vol, all = T)
  siren_vol <- merge(siren_vol,siren_mut_glob_vol, all = T)
  siren_vol <- merge(siren_vol,siren_mut_cnt_vol, all = T)
  siren_vol <- merge(siren_vol,siren_mut_idio_vol, all = T)
  
  siren_vol <- cbind(SIM = rep(i,nrow(siren_vol[VOL != 0])),siren_vol[VOL != 0])
  
  siren_vol <- h_in[siren_vol]
  
  siren_vol_list[[i+1]] <- siren_vol
}



#### Importing firm-level information datasets #####

siren_stat <- read_csv(here::here('data_preparation_sampling','data_summary','data_summary_by_siren.csv'))
siren_stat <- siren_stat[,-1]
setDT(siren_stat)

load(here::here('data_raw','siren_sec.RData'))
siren_sec$SIREN_CHAR <- siren_sec$SIREN
siren_sec$SIREN <- as.numeric(siren_sec$SIREN)

siren_vol_tab0 <- siren_vol_list[[1]]
setDT(siren_vol_tab0)
siren_vol_tab0$SIREN <- as.numeric(siren_vol_tab0$SIREN)

## Joining datasets
setkey(siren_stat,SIREN) # keys setting: SIREN
setkey(siren_vol_tab0,SIREN) # keys setting: SIREN
setkey(siren_sec,SIREN) # keys setting: SIREN

siren_micro_data <- siren_vol_tab0[siren_stat] # joining firm statistics and volatility infos
siren_micro_data <- siren_micro_data%>% filter(SIM ==0) # filtering using only the full dataset
siren_micro_data <- siren_sec[siren_micro_data] # join with information on the main sector of the firm


binning_q <- function(x,q) ## function subsetting a variable into bins
{
  y <- rep(NA,length(x))
  qv <- quantile(x,(1:q/q), na.rm = T)
  qv <- round(qv,4)
  level_vec <- rep(NA,length(qv))
  level_vec[q] <- y[x>qv[q-1]] <- paste0(q,'-quantile - (',qv[q-1],',',qv[q],')')
  level_vec[1] <- y[x<=qv[1]] <- paste0('1-quantile - (',0,',',qv[1],')')
  for(j in 2:(q-1))
  {
    level_vec[j] <- y[x>qv[j-1] & x<=qv[j]] <- paste0(j,'-quantile- (',qv[j-1],',',qv[j],')')
  }
  return(factor(y, levels = level_vec))
}
mod_binning_q <- function(x,q,val) ## function subsetting a variable into bins after excluding specific dummy subsets (e.g binning after excluding all the zeros or ones)
{
  y <- x
  y[x==val] <- 'Dummy'
  y[x!=val] <- binning_q(x[x!=val],q)
  return(y)
}


siren_micro_data <- siren_micro_data %>% mutate(BIN_SIZE = binning_q(as.numeric(AVG_SIZE),100)) %>%
                                          mutate(BIN_HQ = binning_q(as.numeric(HQ),100)) %>%
 #                                         mutate(BIN_DDIQ = binning_q(as.numeric(DDIQ),100)) %>%
                                          mutate(BIN_CQ = mod_binning_q(as.numeric(CQ),100,1)) %>%
                                          mutate(BIN_C0 = binning_q(as.numeric(C0),100)) %>%
                                          mutate(BIN_LOG_SIZE = binning_q(as.numeric(log(AVG_SIZE)),100)) 


######## Producing simple descriptives on diversification (Table 3 and Table 4) ########
library(psych)
library(xtable)
df <- siren_micro_data[,.(MEAN_NUM_Q,S_TOP_MEAN_Q,H_IN_MEAN_Q,CQ)]
descr <- psych::describe(df) 
qs <- c(0.01,0.1,0.90,0.99)
quantiles <- t(apply( df , 2 , quantile , probs = qs , na.rm = TRUE ))
setDT(descr)
table_descr <- cbind(descr[,c('vars','n','mean','sd','skew','kurtosis', 'median')],quantiles)

colnames(table_descr) <- c('',' N Obs. ',' Mean ',' SD ',' Skew. ',' Kurt. ',' Median ',' $1^{st}$ pc. ',' $10^{th}$ pc. ',' $90^{th}$ pc. ',' $99^{th}$ pc.')
rownames(table_descr) <- c('Dest. Mkts.', 'Top Share', 'HH Ind.','Div. Ind.')

fileConn <- file(here::here('results_analysis','tables',"Tab_3_div_summary_statistics.txt"))
writeLines(print(xtable(table_descr, type = "latex")), fileConn)
close(fileConn)



upper_clean <- function(x)
{
  for(i in 1:nrow(x))
  {   
    for(j in i:ncol(x)){x[i,j] <- ''}  
  }
  return(x)
}
cor_df <- siren_micro_data[,.(VOL,MEAN_NUM_Q,S_TOP_MEAN_Q,H_IN_MEAN_Q,CQ)]
cor_mat_q <- upper_clean(format(round(cor(cor_df, use = 'complete.obs'),4), nsmall = 4))
cor_mat_q <- cor_mat_q[!colnames(cor_mat_q)%in%c('VOL'),!colnames(cor_mat_q)%in%c('CQ')]

colnames(cor_mat_q) <-  c('$\\sigma$', 'Dest. Mkts.', 'Top Share', 'HH Ind.')
rownames(cor_mat_q) <-  c('Dest. Mkts.', 'Top Share', 'HH Ind.','Div. Ind.')


fileConn <- file(here::here('results_analysis','tables',"Tab_4_div_correlations.txt"))
writeLines(print(xtable(cor_mat_q, type = "latex")), fileConn)
close(fileConn)


###### Plotting counterfactual distributions (Figure 10)  ###### 
library(ggridges)
library(plyr)

# Call the pdf command to start the plot
pdf(file = here::here('results_analysis','figures','Fig_10_counterfactual_vols.pdf'),   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 8) # The height of the plot in inches

# Create the plot

data_sum_conterfactual <- siren_micro_data[,c('AVG_SIZE',
                                              'MUT_GLOB_VOL',
                                              'MUT_CNT_VOL',
                                              'MUT_IDIO_VOL',
                                              'WGT_VOL_Q')] %>%
                          gather(vol_type,vol_val,2:5) %>% 
                          mutate(BIN_LOG_SIZE = binning_q(as.numeric(log(AVG_SIZE)),5))

p <- data_sum_conterfactual %>%  mutate(vol_type = factor(vol_type,
                                                     levels = c('WGT_VOL_Q','MUT_GLOB_VOL','MUT_CNT_VOL','MUT_IDIO_VOL'))) %>%
                            mutate(vol_type = revalue(vol_type,
                                                      c('MUT_IDIO_VOL' = 'Idio. muted',
                                                        'MUT_GLOB_VOL' = 'Global muted',
                                                        'MUT_CNT_VOL' = 'Dest. muted',
                                                        'WGT_VOL_Q' = 'Total vol.'))) %>%
                            ggplot(aes(x = 0.01 + sqrt(vol_val), y = vol_type)) +
                            scale_x_continuous(limits = c(0.01,2),labels = fmt_decimals2(4)) + 
                            stat_density_ridges(aes(fill = vol_type), alpha = 1.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_grey(labels = c("Dest. specific","Global","Idiosyncratic")) + #,"Volatility on yearly basis"))  +
                            scale_color_grey(labels = c("Dest. specific","Global","Idiosyncratic")) +
                            labs(y = '', x = '', title = '') +  
                            guides(fill=guide_legend(title=""),
                                   col =guide_legend(title="") ) +
                            theme_minimal()  +
                            theme(legend.direction = 'horizontal',
                                  legend.position = 'none',
                                  strip.background = element_blank(),
                                  strip.text.x = element_blank(),
                                  axis.text.x = element_text(size = 20),
                                  axis.text.y = element_text(size = 20)) 
print(p)

dev.off()
detach("package:plyr", unload=TRUE)    
detach("package:ggridges", unload=TRUE)    


### Restrict the data taking the distributions only
restrict_for_summary <- siren_micro_data[,c('MUT_GLOB_VOL','MUT_CNT_VOL','MUT_IDIO_VOL','WGT_VOL_Q','VOL')]
qs_vol <-quantile(restrict_for_summary$WGT_VOL_Q, c(.10,.25,.50,.75,.90), na.rm = T)
qs_glob_mut <- quantile(restrict_for_summary$MUT_GLOB_VOL, c(.10,.25,.50,.75,.90), na.rm = T)
qs_cnt_mut <- quantile(restrict_for_summary$MUT_CNT_VOL, c(.10,.25,.50,.75,.90), na.rm = T)
qs_idio_mut <- quantile(restrict_for_summary$MUT_IDIO_VOL, c(.10,.25,.50,.75,.90), na.rm = T)

sumup <- rbind(qs_vol,
             qs_glob_mut,
             (qs_glob_mut-qs_vol)/qs_vol,
             qs_cnt_mut,
             (qs_cnt_mut-qs_vol)/qs_vol,
             qs_idio_mut,
             (qs_idio_mut-qs_vol)/qs_vol)

fileConn <- file(here::here('results_analysis','tables',"Fig_10_tab_moments_distributions.txt"))
writeLines(print(xtable(sumup, digits=c(4,4,4,4,4,4), type = "latex")), fileConn)
close(fileConn)





#### PLotting distributions diversification (inverse herfindahl) index vs vol #####


### Prepare the dataset
data_sum <- siren_micro_data %>% select(c('BIN_CQ','CQ','VOL','GLOB_VOL','CNT_VOL','IDIO_VOL'))
setDT(data_sum)

data_sum[,CQ := mean(CQ, na.rm = T), keyby = .(BIN_CQ)]

vols <- c('VOL','GLOB_VOL','CNT_VOL','IDIO_VOL')

vol_quant <- 1:9/10 ## define the volatility quantiles to group the firms
data_sum_list <- list()
for(qq in vol_quant)
{
  data_par <- data_sum[ , lapply(.SD, function(x) quantile(x,qq)), keyby = .(BIN_CQ,CQ), .SDcols = vols]
  data_par <- cbind(data_par,BIN_VOL = rep(qq,nrow(data_par)))
  data_sum_list[[which(qq == vol_quant)]] <- data_par
}

data_sum <- do.call(rbind,data_sum_list) %>% gather(vol_type,vol_val,3:6)
setDT(data_sum)


############### Plot of the volatility deciles vs diversification quantiles (Figure 11) ######################

library(plyr)

# Call the pdf command to start the plot
pdf(file = here::here('results_analysis','figures','Fig_11_volatility_vs_div_deciles_all.pdf'),   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 8) # The height of the plot in inches


# Create the plot
p <- data_sum %>%  filter(!is.na(BIN_CQ)) %>% 
              mutate(vol_type = revalue(vol_type, c('CNT_VOL' = 'Dest. specific', 'GLOB_VOL' = 'Common', 'IDIO_VOL'='Idiosyncratic', 'VOL'='Total volatility'))) %>%
              mutate(BIN_CQ = as.numeric(BIN_CQ)) %>% filter(BIN_CQ != 'Dummy') %>% 
              ggplot() +
              geom_point(aes(y = vol_val, x = CQ, group = factor(BIN_VOL),
                             fill = factor(BIN_VOL, levels = 9:1/10), color = factor(BIN_VOL, levels = 9:1/10)),
                         alpha = 0.9, cex = 3) +
              #labels = c("Volatility per Herfindahl index deciles",))  +
              scale_fill_grey() + 
              scale_color_grey() +
              scale_x_continuous(trans = 'log',labels = fmt_decimals2(4)) +
              #labels = c("Volatility per Herfindahl index deciles")) + 
              labs(y = 'Volatility (decile thresholds)', x = 'Diversification (quantile averages - log scale)', title = '') +  
              guides(fill=guide_legend(title=""),
                     col =guide_legend(title="") ) +
              theme_minimal()  +
              theme(legend.direction = 'horizontal',
                    legend.position = 'none',  
                    strip.background = element_blank()) +#,
              #strip.text.x = element_blank()) +
              facet_wrap(~vol_type, ncol =2, scales = 'free')

print(p)
dev.off()

detach("package:plyr", unload=TRUE)    


############### Plot of the volatility deciles vs diversification quantiles restricted ###################
############### to first and last decile of the volatility distribution (i.e. firms with #################
############### most volatile and least volatile growth paths (figure 12) ################################

library(plyr)
p1 <- data_sum %>%  filter(!is.na(BIN_CQ)) %>% filter(BIN_CQ != 'Dummy') %>%  filter(vol_type != 'VOL' ) %>%
  mutate(vol_type = revalue(vol_type, c('CNT_VOL' = 'Dest. specific', 'GLOB_VOL' = 'Common', 'IDIO_VOL'='Idiosyncratic', 'VOL'='Total volatility'))) %>%
  mutate(BIN_CQ = as.numeric(BIN_CQ)) %>% filter( BIN_VOL %in% c(0.9) )%>% filter(vol_type != 'Total volatility' ) %>% 
  ggplot() +
  geom_line(aes(y = vol_val, x = CQ, group = factor(vol_type),
                fill = vol_type, color = vol_type), alpha = 1, cex = 2) +
  #labels = c("Volatility per Herfindahl index deciles",))  +
  scale_fill_grey() + 
  scale_color_grey() +
  scale_x_continuous(trans = 'log',labels = fmt_decimals2(4)) +
  scale_y_continuous(position = 'left') +
  #labels = c("Volatility per Herfindahl index deciles")) + 
  labs(y = 'Volatility (top 10%)', x = 'Diversification (quantile averages - log scale)', title = '') +  
  guides(fill=guide_legend(title=""),
         col =guide_legend(title="") ) +
  theme_minimal()  +
  theme(axis.text.x = element_text(size = 16),
        axis.text.y = element_text(size = 16),
        panel.grid.major.y = element_line(color = 'azure2'),
        panel.grid.major.x = element_blank(),
        panel.background = element_blank(),
        axis.title.x = element_text(size = 18),
        axis.title.y = element_text(size =18),
        legend.direction = 'horizontal',
        legend.position = 'top', 
        legend.title = element_blank(),
        legend.key= element_blank(),
        legend.text = element_text(size = 20),
        strip.text.x = element_text(size = 20, color = "black", face = "bold"),
        strip.text.y = element_text(size = 20, color = "black", face = "bold")) 


p2 <- data_sum %>%  filter(!is.na(BIN_CQ)) %>% filter(BIN_CQ != 'Dummy') %>% #filter(vol_type != 'VOL' ) %>%
  mutate(vol_type = revalue(vol_type, c('CNT_VOL' = 'Dest. specific', 'GLOB_VOL' = 'Common', 'IDIO_VOL'='Idiosyncratic', 'VOL'='Total volatility'))) %>%
  mutate(BIN_CQ = as.numeric(BIN_CQ)) %>% filter( BIN_VOL %in% c(0.1)) %>% filter(vol_type != 'Total volatility' ) %>%
  ggplot() +
  geom_line(aes(y = vol_val, x = CQ, group = factor(vol_type),
                fill = vol_type, color = vol_type), alpha = 1, cex = 2) +
  #labels = c("Volatility per Herfindahl index deciles",))  +
  scale_fill_grey() + 
  scale_color_grey() +
  scale_x_continuous(trans = 'log',labels = fmt_decimals2(4)) +
  scale_y_continuous(position = 'right') +
  #labels = c("Volatility per Herfindahl index deciles")) + 
  labs(y = 'Volatility (bottom 10%)', x = 'Diversification (logs quantile averages)', title = '') +  
  guides(fill=guide_legend(title=""),
         col =guide_legend(title="") ) +
  theme_minimal()  +
  theme(axis.text.x = element_text(size = 16),
        axis.text.y = element_text(size = 16),
        panel.grid.major.y = element_line(color = 'azure2'),
        panel.grid.major.x = element_blank(),
        panel.background = element_blank(),
        axis.title.x = element_text(size = 18),
        axis.title.y = element_text(size =18),
        legend.direction = 'horizontal',
        legend.position = 'top', 
        legend.title = element_blank(),
        legend.key= element_blank(),
        legend.text = element_text(size = 20),
        strip.text.x = element_text(size = 20, color = "black", face = "bold"),
        strip.text.y = element_text(size = 20, color = "black", face = "bold")) 


g_legend<-function(a.gplot){
  tmp <- ggplot_gtable(ggplot_build(a.gplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)}

mylegend<-g_legend(p1)

# Call the pdf command to start the plot
pdf(file = here::here('results_analysis','figures','Fig_12_volatility_vs_div_deciles_top_bottom.pdf'),   # The directory you want to save the file in
    width = 15, # The width of the plot in inches
    height = 8) # The height of the plot in inches


# Create the plot
pp<- grid.arrange(mylegend,arrangeGrob(p1 + theme(legend.position="none"),
                                        p2 + theme(legend.position="none"),
                                        nrow=1),
                   nrow=2,heights=c(1, 10))

print(pp)
dev.off()

detach("package:plyr", unload=TRUE)    


########## Residual regressions and plot versus diversification index (Figure 13) ############
n <- 100

## Prepare the data
siren_micro_data <- siren_micro_data %>%
  mutate(BIN_SIZE = binning_q(as.numeric(AVG_SIZE),n)) %>%
  mutate(BIN_HQ = binning_q(as.numeric(HQ),n)) %>%
  #                mutate(BIN_DDIQ = binning_q(as.numeric(DDIQ),n)) %>%
  mutate(BIN_CQ = mod_binning_q(as.numeric(CQ),n,1)) %>%
  mutate(BIN_C0 = binning_q(as.numeric(C0),n)) %>%
  mutate(BIN_LOG_SIZE = binning_q(as.numeric(log(AVG_SIZE)),n)) 


## Run the regressions
mod_vol0 <- lm(formula = log(VOL) ~ log(AVG_SIZE) + factor(BIN_LOG_SIZE) , data = siren_micro_data)
mod_com0 <- lm(formula = log(GLOB_VOL) ~ log(AVG_SIZE) + factor(BIN_LOG_SIZE) , data = siren_micro_data)
mod_cnt0 <- lm(formula = log(CNT_VOL) ~ log(AVG_SIZE) + factor(BIN_LOG_SIZE) , data = siren_micro_data)
mod_idio0 <- lm(formula = log(IDIO_VOL) ~ log(AVG_SIZE) + factor(BIN_LOG_SIZE) , data = siren_micro_data)


mod_list <- list(list('VOL',mod_vol0),
                 list('GLOB_VOL',mod_com0),
                 list('CNT_VOL',mod_cnt0),
                 list('IDIO_VOL',mod_idio0))

### Store the residuals and collect them into bins
cuts <- 30
data_list <- list()
for(mod in mod_list)
{
  vol_type <- mod[[1]]
  model <- mod[[2]]
  
  output0 <- predict(model,siren_micro_data) %>% 
    cbind(predicted = ., siren_micro_data) %>% 
    mutate(sqr_residual = (log(.data[[vol_type]]) - predicted)^1) %>%
    mutate(res_bin = cut(sqr_residual,cuts)) %>% group_by(res_bin) %>%
    mutate(sqr_res_m = mean(sqr_residual, na.rm = T))
  setDT(output0)
  output0[,CQ_m := mean(CQ, na.rm = T), keyby = .(BIN_CQ)]
  
  
  data_sum_vol <- output0 %>% select(c('BIN_CQ','res_bin','CQ','sqr_res_m')) %>%
    unique() %>% 
    mutate(vol_tp = rep(vol_type, nrow(.)))
  setDT(data_sum_vol)
  
  data_list[[vol_type]] <- data_sum_vol
}

data_sum <- do.call(rbind, data_list)


library(plyr)

# Call the pdf command to start the plot
pdf(file = here::here('results_analysis','figures','Fig_13_residual_plots.pdf'),   # The directory you want to save the file in
    width = 12, # The width of the plot in inches
    height = 8) # The height of the plot in inches


# Create the plot
p <- data_sum %>%  filter(!is.na(BIN_CQ)) %>% 
  filter(BIN_CQ != 'Dummy') %>%
  mutate(vol_tp = revalue(vol_tp, c('CNT_VOL' = 'Dest. specific',
                                    'GLOB_VOL' = 'Common',
                                    'IDIO_VOL'='Idiosyncratic',
                                    'VOL'='Total volatility'))) %>%
  ggplot() +
  geom_point(aes(y = sqr_res_m, x = CQ, color = vol_tp),
             alpha = 0.5, cex = 0.5) +
  #labels = c("Volatility per Herfindahl index deciles",))  +
  scale_fill_grey() + 
  scale_color_grey() +
  scale_x_continuous(trans = 'log',labels = fmt_decimals2(4)) +
  #labels = c("Volatility per Herfindahl index deciles")) + 
  labs(y = 'Regression residuals: ln(Vol.) ~ ln(Size) +  Size Quant.', x = 'Diversification (log scale)', title = '') +  
  guides(fill=guide_legend(title=""),
         col =guide_legend(title="") ) +
  theme_minimal()  +
  theme(axis.text.x = element_text(size = 16),
        axis.text.y = element_text(size = 16),
        panel.grid.major.y = element_line(color = 'azure2'),
        panel.grid.major.x = element_blank(),
        panel.background = element_blank(),
        axis.title.x = element_text(size = 18),
        axis.title.y = element_text(size =18),
        legend.position = 'none', 
        legend.title = element_blank(),
        legend.key= element_blank(),
        legend.text = element_text(size = 20),
        strip.text.x = element_text(size = 20, color = "black", face = "bold"),
        strip.text.y = element_text(size = 20, color = "black", face = "bold"))  +
  facet_wrap(~vol_tp, ncol = 2)

print(p)

dev.off()

detach("package:plyr", unload=TRUE)    


######### Firm-level regressions (not included in the manuscript) ###########
siren_micro_data_decile <- siren_micro_data %>% 
  mutate(BIN_SIZE = binning_q(as.numeric(AVG_SIZE),10)) %>%
  mutate(BIN_HQ = binning_q(as.numeric(HQ),10)) %>%
  # mutate(BIN_DDIQ = binning_q(as.numeric(DDIQ),10)) %>%
  mutate(BIN_CQ = mod_binning_q(as.numeric(CQ),10,1)) %>%
  mutate(BIN_C0 = binning_q(as.numeric(C0),10)) %>%
  mutate(BIN_LOG_SIZE = binning_q(as.numeric(log(AVG_SIZE)),10)) 


mod1 <- lm(formula = log(VOL) ~ CQ, data = siren_micro_data_decile)
mod2 <- lm(formula = log(VOL) ~ CQ + factor(BIN_LOG_SIZE) , data = siren_micro_data_decile)
mod3 <- lm(formula = log(VOL) ~ CQ + factor(BIN_LOG_SIZE) + factor(BIN_LOG_SIZE):factor(BIN_CQ), data = siren_micro_data_decile)
mod4 <- lm(formula = log(VOL) ~ CQ + factor(BIN_LOG_SIZE) + factor(BIN_LOG_SIZE):factor(BIN_CQ) + factor(NC2), data = siren_micro_data_decile)

stargazer(mod1,mod2,mod3,mod4)

mod_com <- lm(formula = log(GLOB_VOL) ~  CQ + factor(BIN_LOG_SIZE) + BIN_LOG_SIZE:BIN_CQ + factor(NC2), data = siren_micro_data_decile)
mod_cnt <- lm(formula = log(CNT_VOL) ~  CQ + factor(BIN_LOG_SIZE) + BIN_LOG_SIZE:BIN_CQ + factor(NC2), data = siren_micro_data_decile)
mod_idio <- lm(formula = log(IDIO_VOL) ~  CQ + factor(BIN_LOG_SIZE) + BIN_LOG_SIZE:BIN_CQ + factor(NC2), data = siren_micro_data_decile)

stargazer(mod_com, mod_cnt, mod_idio)


##### saving main results in R files ######
siren_vol_list0 <-  siren_vol_list[[1]]
save(siren_vol_list0, file=here::here('results_analysis','partial_outputs',"siren_vol_list0.RData"))
