################################################################################
############ Main script for results analysis after model estimation ###########
############ Results for the aggregate volatility - Macro level ################
################################################################################

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

#### Auxiliary packages ####
pacman::p_load(here)
pacman::p_load(data.table)
pacman::p_load(tidyr)
pacman::p_load(readr)
pacman::p_load(dplyr)
pacman::p_load(ggplot2)
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('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'

#### The script reads the data for each sample of firms defined in the data preparation procedure ####

H <- 20 #### number of dataset subsamples

#### 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:H)
{
  cat('\nStarting the step',i,'of the cycle',ifelse(i==0,'(at step 0 the complete dataset is used)\n','\n'))
  
  ######################### Imports and components definition and calculations #########################
  
  cat('\nImporting the outputs of the estimation stage...\n')
  
  
  ### 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)

  ## Import factor loadings
  cat('   Importing the loadings...\n')
  
  load <- setDT(read_csv(here('results_analysis','Inputs',folder_res,paste0('results_',i,'_C.csv')),
                         show_col_types = FALSE))
 
  
  ## 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])
  
  ## 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]
  ## Preliminary operation: identify and store the colunm names of the export value columns
  numeric_cols <- colnames(gr_mat)[3:ncol(gr_mat)]
  ## Force numeric format for columns containing export values
  gr_mat[,(numeric_cols) := lapply(.SD, function(x) as.numeric(x)), .SDcols = numeric_cols]
  
  
  com_comp <- as.matrix(load[,3:ncol(load)])%*%t(as.matrix(fac)) ## common and dest components in one
  
  ## Clean the memory (MIGHT BE NECESSARY)
  ## rm(com_com)
  ## gc()
  
  
  ############################# Factor space reconstruction #################################
  cat('\nFactor space reconstruction..\n')  
  
  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)
  
  fac_list[[i+1]] <- fac_tab

  
  ####################### Aggregate volatility calculations #####################
  
  cat('\nAggregate volatility calculations...\n')
  
  ## Import weights for the time-series selected in the sample (w_sel)
    
  w_sel <- setDT(read_csv(here::here('results_analysis','Inputs',folder_res,folder_samp,paste0('data_sample_w_sel',i,'.csv')),
                          show_col_types = FALSE,
                          name_repair = "minimal"))
  w_sel <- w_sel[,-1]

  
  weights_matrix <- as.matrix(w_sel[,3:(ncol(w_sel)-4)])
  

  ## Rescaling is necessary to refactor the variance and the means of the original time series
  ## that were standardized for the estimation of the factor model
  
  #### Set up the matrix for rescaling
  
  ## Matrix with NA pattern
  w_na <- as.matrix(w_sel[,3:(ncol(w_sel)-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 <- matrix(rep(v_sig,ncol(weights_matrix)), ncol = ncol(weights_matrix))
  mean_mat <- matrix(rep(v_means,ncol(weights_matrix)), ncol = ncol(weights_matrix))
  
  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 
  

  ## Auxiliary functions: sums and means by columns excluding the outliers
  ## NOTICE: if no argument outlier_thresh is passed, outliers are not removed
  colSums_out <- function(x, na.rm, outlier_thresh = Inf){apply(x,2, function(y) sum(y[y>-outlier_thresh&y<outlier_thresh], na.rm=na.rm))}
  colMeans_out <- function(x, na.rm, outlier_thresh = Inf){apply(x,2, function(y) mean(y[y>-outlier_thresh&y<outlier_thresh], na.rm=na.rm))}
  
  ### Aggregations. Three methods of aggregating. 
  ### M1) aggregation through dynamic weights 
  ### M2) taking the means of the series
  ### M3) taking constant weights (at different time points) as in Di Giovanni et al 2012, Econometrica (see the MATLAB function var_frwoth.m in the DiGiovanni_et_al_replication folder)
  ###     This method provides T different aggregates, one for each cross-section (to which a size-based weight is associated)
  
  ### Constructing the aggregate using the actual growth rates (after demeaning)
  
  cat('   Aggregate volatility of the actual growth rates...\n')
  
  gr_actual_mat <- gr_mat[,3:ncol(gr_mat)] - rowMeans(gr_mat[,3:ncol(gr_mat)], na.rm = T)
  
  M1_agg_actual <- colSums_out(weights_matrix*gr_actual_mat, na.rm = T) ## 
  M2_agg_actual_mean <- colMeans_out(gr_actual_mat, na.rm = T)
  m_gr <- as.matrix(gr_actual_mat) 
  M3_matrix_agg_actual <- const_weights_aggregate(weights_matrix,m_gr)
 

  #### Aggregate vectors of each component and their sum (common, destination(cnt) and idiosyncratic)
  
  #### The common global component
  
  ## common global component of the growth rates rescaled including the variance of the original series
  cat('   Aggregate volatility of the global component...\n')
  
  common_global_mat_na <- sd_mat_na*as.matrix(load[,3])%*%t(as.matrix(fac[,1])) 
  common_global_mat <- sd_mat*as.matrix(load[,3])%*%t(as.matrix(fac[,1])) 
  
  ## calculating the three aggregates for the global component
  
  M1_agg_glob <- colSums_out(weights_matrix*common_global_mat, na.rm = T)
  M1_agg_glob_na <- colSums_out(weights_matrix*common_global_mat_na, na.rm = T)
  M2_agg_glob_mean <- colMeans_out(common_global_mat, na.rm = T)
  M2_agg_glob_mean_na <- colMeans_out(common_global_mat_na, na.rm = T)
  M3_matrix_agg_glob <- const_weights_aggregate(weights_matrix,common_global_mat)
  M3_matrix_agg_glob_na <- const_weights_aggregate(weights_matrix,common_global_mat_na)
  

  ## common destination component of the growth rates rescaled including the variance of the original series
  cat('   Aggregate volatility of the destination component...\n')
  
  common_destination_mat_na <- sd_mat_na*as.matrix(load[,4:ncol(load)])%*%t(as.matrix(fac[,2:ncol(fac)])) 
  common_destination_mat <-  sd_mat*as.matrix(load[,4:ncol(load)])%*%t(as.matrix(fac[,2:ncol(fac)])) 
  
  
  ## calculating the three aggregates for the destination-specific component
  M1_agg_dest <- colSums_out(weights_matrix*common_destination_mat, na.rm = T)
  M1_agg_dest_na <- colSums_out(weights_matrix*common_destination_mat_na, na.rm = T)
  M2_agg_dest_mean <- colMeans_out(common_destination_mat, na.rm = T)
  M2_agg_dest_mean_na <- colMeans_out(common_destination_mat_na, na.rm = T)
  M3_matrix_agg_dest <- const_weights_aggregate(weights_matrix,common_destination_mat)
  M3_matrix_agg_dest_na <- const_weights_aggregate(weights_matrix,common_destination_mat_na)
  
  
  ## sum of the global and dest component of the growth rates rescaled including the variance of the original series
  
  cat('   Aggregate volatility of the common component...\n')
  
  common_mat_na <-  sd_mat_na*as.matrix(load[,3:ncol(load)])%*%t(as.matrix(fac[,1:ncol(fac)]))
  common_mat <-  sd_mat*as.matrix(load[,3:ncol(load)])%*%t(as.matrix(fac[,1:ncol(fac)]))
  
  ## calculating the three aggregates for the destination-specific component
  M1_agg_common <- colSums_out(weights_matrix*common_mat, na.rm = T)
  M1_agg_common_na <- colSums_out(weights_matrix*common_mat_na, na.rm = T)
  M2_agg_common_mean <- colMeans_out(common_mat, na.rm = T)
  M2_agg_common_mean_na <- colMeans_out(common_mat_na, na.rm = T)
  M3_matrix_agg_common <- const_weights_aggregate(weights_matrix,common_mat)
  M3_matrix_agg_common_na <- const_weights_aggregate(weights_matrix,common_mat_na)
  
  
  ## idiosyncratic component of the growth rates calculated by difference from the actual growth rates and rescaled common component
  cat('Aggregate volatility of the idiosyncratic component...\n')
  
  
  idio_mat <- as.matrix(gr_actual_mat - common_mat) 
  ## calculating the three aggregates for the idio component
  M1_agg_idio <- colSums_out(weights_matrix*idio_mat, na.rm = T) 
  M2_agg_idio_mean <- colMeans_out(idio_mat, na.rm = T)
  M3_matrix_agg_idio <- const_weights_aggregate(weights_matrix,idio_mat)
  
  

  #### Relative volatility calculations
  
  M1_agg_components_matrix <- cbind(M1_agg_common, M1_agg_common_na,
                                    M1_agg_glob, M1_agg_glob_na,
                                    M1_agg_dest, M1_agg_dest_na,
                                    M1_agg_idio)
  
  M1_corr <- as.data.frame(cbind(cor_comm = cor(M1_agg_actual, M1_agg_common),
                                 cor_glob = cor(M1_agg_actual, M1_agg_glob), 
                                 cor_dest = cor(M1_agg_actual ,M1_agg_dest), 
                                 cor_idio = cor(M1_agg_actual, M1_agg_idio)))
  
  M2_mean_components_matrix <- cbind(M2_agg_common_mean, M2_agg_common_mean_na,
                                     M2_agg_glob_mean, M2_agg_glob_mean_na,
                                     M2_agg_dest_mean, M2_agg_dest_mean_na,
                                     M2_agg_idio_mean)
  
  M2_corr <- as.data.frame(cbind(cor_comm = cor(M2_agg_actual_mean, M2_agg_common_mean),
                                 cor_glob = cor(M2_agg_actual_mean, M2_agg_glob_mean),
                                 cor_dest = cor(M2_agg_actual_mean, M2_agg_dest_mean),
                                 cor_idio = cor(M2_agg_actual_mean, M2_agg_idio_mean)))
  
  
  M3_components_matrix <- cbind(M3_agg_common_prime = colMeans(M3_matrix_agg_common$var_mat),
                                M3_agg_common_na_prime = colMeans(M3_matrix_agg_common_na$var_mat),
                                M3_agg_glob_prime = colMeans(M3_matrix_agg_glob$var_mat), 
                                M3_agg_glob_na_prime = colMeans(M3_matrix_agg_glob_na$var_mat), 
                                M3_agg_dest_prime = colMeans(M3_matrix_agg_dest$var_mat),
                                M3_agg_dest_na_prime = colMeans(M3_matrix_agg_dest_na$var_mat),
                                M3_agg_idio_prime = colMeans(M3_matrix_agg_idio$var_mat))
  
  M3_corr <- as.data.frame(cbind(cor_comm = cor(colMeans(M3_matrix_agg_actual$var_mat), colMeans(M3_matrix_agg_common$var_mat)),
                                 cor_glob = cor(colMeans(M3_matrix_agg_actual$var_mat), colMeans(M3_matrix_agg_glob$var_mat)),
                                 cor_dest = cor(colMeans(M3_matrix_agg_actual$var_mat), colMeans(M3_matrix_agg_dest$var_mat)),
                                 cor_idio = cor(colMeans(M3_matrix_agg_actual$var_mat), colMeans(M3_matrix_agg_idio$var_mat))))
  

  
  cat('\nStoring the ouputs...\n')
  
  
  agg_cor_list[[i+1]] <- rbind(M1_correlations = M1_corr, M2_correlations = M2_corr, M3_correlations = M3_corr)
  
  
  if(i == 0)
  {
   
    dates <- as.character(seq(as.Date("1994/1/1"), as.Date("2017/10/1"), by = "quarter"))
    
    recap_dataset <- rbind(cbind(source = 'DFM Complete Dataset', type = 'Volatility',  # common vol time series
                                 component = 'Actual', period = dates, 
                                 value = M3_matrix_agg_actual$sd, 
                                 std_err = as.vector(M3_matrix_agg_actual$se_sd)),
                           cbind(source = 'DFM Complete Dataset', type = 'Volatility',  # common vol time series
                                 component = 'Common', period = dates, 
                                 value = M3_matrix_agg_common_na$sd, 
                                 std_err = as.vector(M3_matrix_agg_common_na$se_sd)),
                           cbind(source = 'DFM Complete Dataset', type = 'Volatility',  # global vol time series
                                 component = 'Global', period = dates, 
                                 value = M3_matrix_agg_glob_na$sd, 
                                 std_err = as.vector(M3_matrix_agg_glob_na$se_sd)),
                           cbind(source = 'DFM Complete Dataset', type = 'Volatility', # destination vol time series
                                 component = 'Destination', period = dates, 
                                 value = M3_matrix_agg_dest_na$sd, 
                                 std_err = as.vector(M3_matrix_agg_dest_na$se_sd)),
                           cbind(source = 'DFM Complete Dataset', type = 'Volatility', # idio vol time series
                                 component = 'Idio', period = dates, 
                                 value = M3_matrix_agg_idio$sd, 
                                 std_err = as.vector(M3_matrix_agg_idio$se_sd))) 
    
    
    
    recap_dataset <- as.data.frame(recap_dataset)
    rownames(recap_dataset) <- NULL
    vol_comparison_ts_list[['full_dataset']] <- recap_dataset  
    
    matrix_list[['full_dataset']][['common']] <- common_mat_na
    matrix_list[['full_dataset']][['idio']] <- idio_mat
    
    ###### Calculating the bootstrap standard error associated to the volatility point estimates on the whole time span
    
    bs_std_err <- function(mat,n,label) ## bs standard error of the volatility ratio for the matrix mat
    {  
      bs_std_err_95 <- bootstrap_CI(sapply(1:n,
                                           function(x) mean(sqrt(diag(mat[,sample(1:ncol(mat))])),na.rm = T))/mean(M3_matrix_agg_actual$sd),
                                    0.95)
      bs_std_err_05 <- bootstrap_CI(sapply(1:n,
                                           function(x) mean(sqrt(diag(mat[,sample(1:ncol(mat))])),na.rm = T))/mean(M3_matrix_agg_actual$sd),
                                    0.05)
      bs_std_err_99 <- bootstrap_CI(sapply(1:n,
                                           function(x) mean(sqrt(diag(mat[,sample(1:ncol(mat))])),na.rm = T))/mean(M3_matrix_agg_actual$sd),
                                    0.99)
      bs_std_err_01 <- bootstrap_CI(sapply(1:n,
                                           function(x) mean(sqrt(diag(mat[,sample(1:ncol(mat))])),na.rm = T))/mean(M3_matrix_agg_actual$sd),
                                    0.01)
      out <- c(bs_se_01 = bs_std_err_01, bs_se_05 = bs_std_err_05, bs_se_95 = bs_std_err_95,bs_se_99 = bs_std_err_99)
      names(out) <- paste0('M3_',label,'_',names(out))
      return(t(as.data.frame(out)))
    }

    
    agg_results <- cbind(SIM = c(as.integer(i)),
                         sqrt(t(diag(cov(M1_agg_components_matrix)/var(M1_agg_actual)))),
                         sqrt(t(diag(cov(M2_mean_components_matrix)/var(M2_agg_actual_mean)))),
                         sqrt(t(diag(cov(M3_components_matrix)/var(colMeans(M3_matrix_agg_actual$var_mat))))),
                         M3_agg_common = mean(M3_matrix_agg_common$sd)/mean(M3_matrix_agg_actual$sd),
                         bs_std_err(M3_matrix_agg_common$var_mat,10000,'common'),
                         M3_agg_common_na = mean(M3_matrix_agg_common_na$sd)/mean(M3_matrix_agg_actual$sd),
                         bs_std_err(M3_matrix_agg_common_na$var_mat,10000,'common_na'),
                         M3_agg_glob = mean(M3_matrix_agg_glob$sd)/mean(M3_matrix_agg_actual$sd),
                         bs_std_err(M3_matrix_agg_glob$var_mat,10000,'global'),
                         M3_agg_glob_na = mean(M3_matrix_agg_glob_na$sd)/mean(M3_matrix_agg_actual$sd),
                         bs_std_err(M3_matrix_agg_glob_na$var_mat,10000,'global_na'),
                         M3_agg_dest = mean(M3_matrix_agg_dest$sd)/mean(M3_matrix_agg_actual$sd),
                         bs_std_err(M3_matrix_agg_dest$var_mat,10000,'dest'),
                         M3_agg_dest_na = mean(M3_matrix_agg_dest_na$sd)/mean(M3_matrix_agg_actual$sd),
                         bs_std_err(M3_matrix_agg_dest$var_mat,10000,'dest_na'),
                         M3_agg_idio = mean(M3_matrix_agg_idio$sd)/mean(M3_matrix_agg_actual$sd),
                         bs_std_err(M3_matrix_agg_idio$var_mat,10000,'idio'))
    
    
    agg_vol_list[[i+1]] <- agg_results


  }else{
    
  agg_results <- cbind(SIM = c(as.integer(i)),
                       sqrt(t(diag(cov(M1_agg_components_matrix)/var(M1_agg_actual)))),
                       sqrt(t(diag(cov(M2_mean_components_matrix)/var(M2_agg_actual_mean)))),
                       sqrt(t(diag(cov(M3_components_matrix)/var(colMeans(M3_matrix_agg_actual$var_mat))))),
                       M3_agg_common = mean(M3_matrix_agg_common$sd)/mean(M3_matrix_agg_actual$sd),
                       M3_agg_common_na = mean(M3_matrix_agg_common_na$sd)/mean(M3_matrix_agg_actual$sd),
                       M3_agg_glob = mean(M3_matrix_agg_glob$sd)/mean(M3_matrix_agg_actual$sd),
                       M3_agg_glob_na = mean(M3_matrix_agg_glob_na$sd)/mean(M3_matrix_agg_actual$sd),
                       M3_agg_dest = mean(M3_matrix_agg_dest$sd)/mean(M3_matrix_agg_actual$sd),
                       M3_agg_dest_na = mean(M3_matrix_agg_dest_na$sd)/mean(M3_matrix_agg_actual$sd),
                       M3_agg_idio = mean(M3_matrix_agg_idio$sd)/mean(M3_matrix_agg_actual$sd))
  
    
    
    agg_vol_list[[i+1]] <- agg_results
    
    
  dates <- as.character(seq(as.Date("1994/1/1"), as.Date("2017/10/1"), by = "quarter"))
    
  recap_dataset <- rbind(cbind(source = 'DFM Sample Datasets', type = 'Volatility',  # common vol time series
                                 component = 'Actual', period = dates, 
                                 value = M3_matrix_agg_actual$sd),
                           cbind(source = 'DFM Sample Datasets', type = 'Volatility',  # common vol time series
                                 component = 'Common', period = dates, 
                                 value = M3_matrix_agg_common_na$sd),
                           cbind(source = 'DFM Sample Datasets', type = 'Volatility',  # global vol time series
                                 component = 'Global', period = dates, 
                                 value = M3_matrix_agg_glob_na$sd),
                           cbind(source = 'DFM Sample Datasets', type = 'Volatility', # destination vol time series
                                 component = 'Destination', period = dates, 
                                 value = M3_matrix_agg_dest_na$sd),
                           cbind(source = 'DFM Sample Datasets', type = 'Volatility', # idio vol time series
                                 component = 'Idio', period = dates, 
                                 value = M3_matrix_agg_idio$sd)) 
    
    recap_dataset <- as.data.frame(cbind(SIM= i, recap_dataset))
    rownames(recap_dataset) <- NULL
    vol_sample_ts_list[[i]] <- recap_dataset
  }
  
  
  
  ######################## Destination-level volatility  calculations ###################
  
  cat('\nDestination-level volatility calculations...\n')
  
  
  cnts <- colnames(load)[4:ncol(load)] ## save the country names in a vector
  cnt_comps <- list()
  for(s in cnts)
  {
    ind <- which(load[,s,with = F]!=0)   ## store the indeces in the loading matrix 
    ## to restrict the calculation to the local series
    
    cnt_comps[[s]] <- list(as.matrix(load[ind,c('common',s), with = F])%*%t(as.matrix(fac[,c('common',s), with = F])),
                           as.matrix(load[ind,c(s), with = F])%*%t(as.matrix(fac[,c(s), with = F])),
                           ind)
  }  
  
  
  agg_vol_fun <- function(x) ## function for the aggregation at the destination level
  {
    ## The function is designed to be applied to the elements of the list cnt_comps
    ## x is an element of the list cnt_comps
    ## x[[1]] is the matrix of common+destination component restricted to the destination series
    ## x[[2]] is the matrix of destination components restricted to the destination series
    ## x[[3]] is a vector with the information of the index of the local series for each destination
    
    min <- min(x[[3]]) 
    max <- max(x[[3]])
    n_f <- length(x[[3]])
    if(n_f==1) ## if there is only one time-series (this case should be excluded in the data preparation and in the estimation step)
    {return(c(n_f,rep(NA,6)))}else{  
      
      ## Matrix with NA pattern
      w_na <- as.matrix(w_sel[x[[3]],3:(ncol(w_sel)-4)])
      w_na[!is.na(w_na)] <- 1
      
      restr_weights_matrix <- weights_matrix[x[[3]],]
      
      denom_weights_matrix <- matrix(rep(rowSums(t(restr_weights_matrix), na.rm = T),
                                         ncol(t(restr_weights_matrix))),
                                     ncol = ncol(t(restr_weights_matrix)))
      
      dest_weights_matrix <- restr_weights_matrix/t(denom_weights_matrix)
      
      ## Matrix with time-series specific standard deviations and means
      ## (the NA pattern of the original matrix is imposed already at this stage)
      dest_sd_mat_na <- matrix(rep(v_sig[x[[3]]],ncol(dest_weights_matrix)), ncol = ncol(dest_weights_matrix))*w_na  
      dest_mean_mat_na <- matrix(rep(v_means[x[[3]]],ncol(dest_weights_matrix)), ncol = ncol(dest_weights_matrix))*w_na 
      
      
      ## common destination component of the growth rates rescaled including the variance of the original series
      common_destination_mat_na <-  dest_sd_mat_na*x[[1]]
      
      
      ## calculating the three aggregates for the destination-specific component
      M1_agg_com_dest_na <- colSums_out(dest_weights_matrix*common_destination_mat_na, na.rm = T)
      M2_agg_com_dest_mean_na <- colMeans_out(common_destination_mat_na, na.rm = T)
      M3_matrix_agg_com_dest_na <- const_weights_aggregate(dest_weights_matrix,common_destination_mat_na)$sd
      
      
      ## sum of the global and dest component of the growth rates rescaled including the variance of the original series
      destination_mat_na <-  dest_sd_mat_na*x[[2]]
      
      ## calculating the three aggregates for the destination-specific component
      M1_agg_dest_na <- colSums_out(dest_weights_matrix*destination_mat_na, na.rm = T)
      M2_agg_dest_mean_na <- colMeans_out(destination_mat_na, na.rm = T)
      M3_matrix_agg_dest_na <- const_weights_aggregate(dest_weights_matrix,destination_mat_na)$sd
      
      
      return(c(n_f,
               sd(M1_agg_dest_na, na.rm = T),
               sd(M2_agg_dest_mean_na, na.rm = T),
               mean(M3_matrix_agg_dest_na),
               sd(M1_agg_com_dest_na, na.rm = T),
               sd(M2_agg_com_dest_mean_na, na.rm = T),
               mean(M3_matrix_agg_com_dest_na)))
      
    }  
  }
  
  
  cnt_tab <- as.data.frame(cbind(names(cnt_comps),
                                 do.call(rbind,lapply(cnt_comps,agg_vol_fun))))
  
  colnames(cnt_tab) <- c('PYOD', 'FLOWS', 'VOL_DYN_WGT', 'VOL_MEANS',
                         'VOL_CONST_WGT','VOL_DYN_WGT_RES', 'VOL_MEANS_RES',
                         'VOL_CONST_WGT_RES')
  dest_vol_list[[i+1]] <- cnt_tab

}

agg_vol_tab_samples <- as.data.frame(do.call(rbind,agg_vol_list[-1]))
agg_vol_tab_complete_dataset <- as.data.frame(do.call(rbind,agg_vol_list[1]))

setDT(agg_vol_tab_samples)
setDT(agg_vol_tab_complete_dataset)

####################### Preparing the comparison with the volatility estimates from Di Giovanni and coauthors ####################

replication_vol_list <- list()
replication_cor_list <- list()


models <- c('quarterly','yearly','destsec_quarterly','destsec_yearly') ### List of the models that are replicated
for(model in models)
{
  cat('Replication of Di Giovanni et al. (2014) on our dataset. Model reference:',model ,'\n')
  
    ##### Import the dataset
    rep_dataset  <- setDT(read_csv(here('DiGiovanni_et_al_replication',
                                        'Output',
                                        paste0('growth_regs_capays_',model,'.csv')),
                                   show_col_types = FALSE))
    colnames(rep_dataset)[which(colnames(rep_dataset)=='year')] <- 'period_numeric'
    colnames(rep_dataset)[which(colnames(rep_dataset)=='period')] <- 'period_numeric'
    
    #### Transform the dataset from the long to the wide form for any of the variables of interest: 
    #### i) matrix of weights, ii) actual growth rates, iii) idio component, iv) macro component
    
    convert_period <- function(x){paste0(floor(x/100),'.0',((x/100-floor(x/100))*4)+1)}
    rep_dataset[,period_numeric := convert_period(period_numeric)]
    if(grepl('destsec',model))
    {
    rep_dataset[,gr_macro := gr_macro + gr_sec]  
    }  
    
    rep_weights_matrix <- dcast(rep_dataset,id_fp ~ period_numeric,fun.agg  = function(x) sum(x,na.rm=T), value.var = "w_fp")
    rep_gr_matrix <- dcast(rep_dataset,id_fp ~ period_numeric,fun.agg  = function(x) sum(x,na.rm=T), value.var = "gr_capays")
    rep_idio_matrix <- dcast(rep_dataset,id_fp ~ period_numeric,fun.agg  = function(x) sum(x,na.rm=T), value.var = "gr_idio")
    rep_com_matrix <- dcast(rep_dataset,id_fp ~ period_numeric,fun.agg  = function(x) sum(x,na.rm=T), value.var = "gr_macro")

    
    f_nan(rep_weights_matrix)
    f_nan(rep_gr_matrix)
    f_nan(rep_idio_matrix)
    f_nan(rep_com_matrix)
    
    selection_flows <- rowSums(is.na(rep_gr_matrix))<= ncol(rep_gr_matrix)-4
    
    rep_weights_matrix <- rep_weights_matrix[selection_flows]
    rep_gr_matrix <- rep_gr_matrix[selection_flows]
    rep_idio_matrix <- rep_idio_matrix[selection_flows]
    rep_com_matrix <- rep_com_matrix[selection_flows]
    
    weights_matrix <- as.matrix(rep_weights_matrix[,-1])
    
    
    ### Constructing the aggregate using the actual growth rates (after demeaning)
    gr_actual_mat <- as.matrix(rep_gr_matrix[,-1]) - rowMeans(rep_gr_matrix[,-1], na.rm = T)
    
    M1_agg_actual <- colSums_out(weights_matrix*gr_actual_mat, na.rm = T)
    M2_agg_actual_mean <- colMeans_out(gr_actual_mat, na.rm = T)
    M3_matrix_agg_actual <- const_weights_aggregate(weights_matrix,gr_actual_mat)
    
    #### Aggregate vectors of each component and their sum (common, destination(cnt) and idiosyncratic)
    
    #### The common global component
    
    ## common global component of the growth rates rescaled including the variance of the original series
    common_mat <-  as.matrix(rep_com_matrix[,-1])
    
    ## calculating the three aggregates for the global component
    M1_agg_common <- colSums_out(weights_matrix*common_mat, na.rm = T)
    M2_agg_common_mean <- colMeans_out(common_mat, na.rm = T)
    M3_matrix_agg_com <- const_weights_aggregate(weights_matrix,common_mat)
    
    ## idiosyncratic component of the growth rates calculated by difference from the actual growth rates and rescaled common component
    idio_mat <- as.matrix(rep_idio_matrix[,-1]) 
    ## calculating the three aggregates for the idio component
    M1_agg_idio <- colSums_out(weights_matrix*idio_mat, na.rm = T) 
    M2_agg_idio_mean <- colMeans_out(idio_mat, na.rm = T)
    M3_matrix_agg_idio <- const_weights_aggregate(weights_matrix,idio_mat)
    
    
    ## Construct the sum of the aggregate components 
    agg_comp_mean <- colMeans_out(common_mat + idio_mat, na.rm = T)
    
    #### Aggregate relative volatility calculations
    
    M1_agg_components_matrix <- cbind(M1_agg_common, M1_agg_idio)
    M1_corr <- as.data.frame(cbind(cor_com = cor(M1_agg_actual, M1_agg_common), 
                                   cor_idio = cor(M1_agg_actual, M1_agg_idio)))
    
    M2_mean_components_matrix <- cbind(M2_agg_common_mean, M2_agg_idio_mean)
    M2_corr <- as.data.frame(cbind(cor_com = cor(M2_agg_actual_mean, M2_agg_common_mean),
                                   cor_idio = cor(M2_agg_actual_mean, M2_agg_idio_mean)))
    
    M3_components_matrix <- cbind(M3_agg_common_prime = colMeans(M3_matrix_agg_com$var_mat),
                                  M3_agg_idio_prime = colMeans(M3_matrix_agg_idio$var_mat))
    
    M3_corr <- as.data.frame(cbind(cor_com = cor(colMeans(M3_matrix_agg_actual$var_mat), colMeans(M3_matrix_agg_com$var_mat)),
                                   cor_idio = cor(colMeans(M3_matrix_agg_actual$var_mat), colMeans(M3_matrix_agg_idio$var_mat))))
    
    bs_std_err <- function(mat,n,label) ## bs standard error of the volatility ratio for the matrix mat
    {  
      bs_std_err_95 <- bootstrap_CI(sapply(1:n,
                                           function(x) mean(sqrt(diag(mat[,sample(1:ncol(mat))])),na.rm = T))/mean(M3_matrix_agg_actual$sd),
                                    0.95)
      bs_std_err_05 <- bootstrap_CI(sapply(1:n,
                                           function(x) mean(sqrt(diag(mat[,sample(1:ncol(mat))])),na.rm = T))/mean(M3_matrix_agg_actual$sd),
                                    0.05)
      bs_std_err_99 <- bootstrap_CI(sapply(1:n,
                                           function(x) mean(sqrt(diag(mat[,sample(1:ncol(mat))])),na.rm = T))/mean(M3_matrix_agg_actual$sd),
                                    0.99)
      bs_std_err_01 <- bootstrap_CI(sapply(1:n,
                                           function(x) mean(sqrt(diag(mat[,sample(1:ncol(mat))])),na.rm = T))/mean(M3_matrix_agg_actual$sd),
                                    0.01)
      out <- c(bs_se_01 = bs_std_err_01, bs_se_05 = bs_std_err_05, bs_se_95 = bs_std_err_95,bs_se_99 = bs_std_err_99)
      names(out) <- paste0('M3_',label,'_',names(out))
      return(t(as.data.frame(out)))
    }
    replication_results <- cbind(SIM = paste('Replication -',model),
                                 sqrt(t(diag(cov(M1_agg_components_matrix)/var(M1_agg_actual)))),
                                 sqrt(t(diag(cov(M2_mean_components_matrix)/var(M2_agg_actual_mean)))),
                                 sqrt(t(diag(cov(M3_components_matrix)/var(colMeans(M3_matrix_agg_actual$var_mat))))),
                                 M3_agg_common = mean(M3_matrix_agg_com$sd)/mean(M3_matrix_agg_actual$sd),
                                 bs_std_err(M3_matrix_agg_com$var_mat,10000,'com'),
                                 M3_agg_idio = mean(M3_matrix_agg_idio$sd)/mean(M3_matrix_agg_actual$sd),
                                 bs_std_err(M3_matrix_agg_com$var_mat,10000,'idio'))
    
    
  
    
    replication_vol_list[[paste0('Replication - ',model)]] <- as.data.frame(replication_results)
    replication_cor_list[[paste0('Replication - ',model)]] <- as.data.frame(rbind(M1_corr = M1_corr,
                                                                                  M2_corr = M2_corr,
                                                                                  M3_corr = M3_corr))
    if(grepl('yearly',model))
    {
    dates <- as.character(seq(as.Date("1994/1/1"), as.Date("2017/1/1"), by = "year"))
    }else{
    dates <- as.character(seq(as.Date("1994/1/1"), as.Date("2017/10/1"), by = "quarter"))  
    } 
    
    recap_dataset <- rbind(cbind(source = paste('Replication -', model), type = 'Volatility',  # common vol time series
                                 component = 'Actual', period = dates, 
                                 value = M3_matrix_agg_actual$sd, 
                                 std_err = as.vector(M3_matrix_agg_actual$se_sd)),
                           cbind(source = paste('Replication -', model), type = 'Volatility',  # common vol time series
                                 component = 'Common', period = dates, 
                                 value = M3_matrix_agg_com$sd, 
                                 std_err = as.vector(M3_matrix_agg_com$se_sd)),
                           cbind(source = paste('Replication -', model), type = 'Volatility', # idio vol time series
                                 component = 'Idio', period = dates, 
                                 value = M3_matrix_agg_idio$sd, 
                                 std_err = as.vector(M3_matrix_agg_idio$se_sd))) 
    
    recap_dataset <- as.data.frame(recap_dataset)
    rownames(recap_dataset) <- NULL
    vol_comparison_ts_list[[paste0('rep_',model)]] <- recap_dataset  
    
    matrix_list[[paste0('rep_',model)]][['common']] <- common_mat
    matrix_list[[paste0('rep_',model)]][['idio']] <- idio_mat
    
}

replication_vol_tab <- do.call(rbind, replication_vol_list)
setDT(replication_vol_tab)


###### Joining the results on destination-specific volatility ######
cnt_vols_mod <- list()
for(i in 2:(H+1))
{
  cnt_vols_mod[[i]] <- cbind(SIM=rep(i-1,nrow(dest_vol_list[[i]])),dest_vol_list[[i]])  
}

cnt_vols_tab <- do.call(rbind,cnt_vols_mod)
cnt_vols_tab$FLOWS <- as.integer(cnt_vols_tab$FLOWS)



#########################################################
########### Plots and tables for Section 4.1 ############
#########################################################

##### Plot: common factor and FRED comparison ######

factor_table <- do.call(rbind,fac_list)  ## transform the factors list in a table

french_export <- read_csv(here('results_analysis',"french_export.csv")) ## import the french export

french_export$Value <- french_export$Value - mean(french_export$Value)


identification <- function(x,y){if(sum((x-y)^2)<sum((-x-y)^2)){x}else{-x}}
identification1 <- function(x){if(x[62]<0){x}else{-x}}

factor_pval <- factor_table %>% group_by(SIM,FACTOR) %>%
                                mutate(VALUE = identification1(VALUE)) %>% ungroup() %>%
                                group_by(TIME,FACTOR) %>% 
                                mutate(AVG = mean(VALUE))  %>%
                                mutate(bts_cmax90 = bootstrap_CI(VALUE,0.95)) %>%
                                mutate(bts_cmin90 = bootstrap_CI(VALUE,0.05)) %>%
                                mutate(bts_cmax99 = bootstrap_CI(VALUE,0.99)) %>%
                                mutate(bts_cmin99 = bootstrap_CI(VALUE,0.01)) %>%
                                mutate(MIN =  min(VALUE)) %>%
                                mutate(MAX = max(VALUE)) %>%
                                mutate(CI=ci_fun(VALUE))%>%
                                mutate(p_CI = CI/AVG) %>%
                                ungroup() 


french_export$Date <- as.Date(french_export$Date,"%d/%m/%Y")

quarters <- as.data.frame(cbind(TIME = 1:96,Quarter = french_export$Date))

setDT(factor_pval)
setDT(quarters)

setkey(quarters,TIME)
setkey(factor_pval,TIME)

factor_pval <- quarters[factor_pval]

min_d <- as.Date("1994-1","%Y-%m")
max_d <- as.Date("2017-10","%Y-%m")
factor_pval$Quarter <- as.Date(factor_pval$Quarter, origin = "1970-01-01")


####### Factor plot with factor and FRED (Figure 5)  #######


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

# Create the plot
p<-factor_pval[,c('Quarter','FACTOR','AVG','bts_cmin99','bts_cmax99','bts_cmin90','bts_cmax90')] %>% 
            unique  %>% filter(FACTOR %in% c('common'))%>% cbind(french_export) %>%
            ggplot() +
            geom_ribbon(aes(x = Quarter, ymin = 3*(AVG + bts_cmin99), ymax = 3*(AVG + bts_cmax99), group = FACTOR,fill = "blue"),
                        alpha = 0.3)+
            geom_ribbon(aes(x = Quarter, ymin = 3*(AVG + bts_cmin90), ymax = 3*(AVG + bts_cmax90), group = FACTOR,fill = "red"),
                        alpha = 0.3)+  
            geom_line(aes(x = Quarter, y = 3*AVG, group = FACTOR,color = 'black', linetype = 'dashed')) +
            theme_light() +
            theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 12),
                  axis.text.y = element_text(size=12),
                  panel.grid.minor.x = element_blank())+
            #panel.grid.minor.y = element_blank())+
            labs(y = '', x = '', title = '') +
            geom_line(aes(x = Quarter, y = Value, color = 'darkred', linetype = 'solid')) +
            scale_x_date(limits=c(min_d,max_d), date_labels = "%b/%Y", breaks = scales::pretty_breaks(n = 9)) + 
            scale_y_continuous(breaks = scales::pretty_breaks(n = 5)) + 
            scale_fill_manual("Confidence intervals", labels = c("99% CI (bootstrap)","90% CI (bootstrap)"),   values=c("blue"="blue","red"="red"))+ 
            scale_color_manual("Series", labels = c("Global factor","FRED agg."), values=c("black"="black","darkred"="darkred")) +
            scale_linetype_manual("Series", labels = c("Global factor","FRED agg."), values=c("solid"="dashed","dashed"="solid"))

print(p)
# Run dev.off() to create the file!
dev.off()


###### Selected destination factors #######

pdf(file = here('results_analysis','figures','Sel_dest_factors.pdf'),   # The directory you want to save the file in
    width = 18, # The width of the plot in inches
    height = 8) # The height of the plot in inches


p <- factor_pval[,c('Quarter','FACTOR','AVG','bts_cmin99','bts_cmax99','bts_cmin90','bts_cmax90')] %>% 
  unique   %>%  filter(FACTOR %in% c('IT','AT','CN','US','EE','LT')) %>% 
  ggplot() +
  geom_ribbon(aes(x = Quarter, ymin = AVG + bts_cmin99, ymax = AVG + bts_cmax99, group = FACTOR,fill = "blue"), alpha = 0.2)+
  geom_ribbon(aes(x = Quarter, ymin = AVG + bts_cmin90, ymax = AVG + bts_cmax90, group = FACTOR,fill = "red"),alpha = 0.2)+  
  geom_line(aes(x = Quarter, y = AVG, group = FACTOR)) +
  # geom_ribbon(aes(x = Time, ymin = min, ymax = max, group = Factor, fill = "black"), alpha = 0.1)+
  theme_light() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  labs(y = '', x = '', title = '') +
  scale_x_date(limits=c(min_d,max_d), date_labels = "%b/%Y", breaks = scales::pretty_breaks(n = 10)) +
  scale_fill_manual("Confidence intervals", labels = c("99% CI (bootstrap)","90% CI (bootstrap)"),   values=c("blue"="blue","red"="red")) +
  facet_wrap(~FACTOR, ncol = 3)
print(p)

# Run dev.off() to create the file!
dev.off()


##### Plot: DFMs vs SODMs time series volatility comparison (Figure 7) ####
    
min_d <- as.Date("1994-1-1","%Y-%m-%d")
max_d <- as.Date("2017-10-31","%Y-%m-%d")

vol_comparison_ts_tab <- do.call(rbind,vol_comparison_ts_list)

volatility_comparison_plot_tab <- vol_comparison_ts_tab %>% 
                                  filter((component %in% c('Common', 'Idio'))) %>% 
                                  filter(source != 'Replication - destsec_yearly') 

names_labels <- unique(volatility_comparison_plot_tab$source)
source_labels <- c("DFM - Global and destination factors - Quarterly", "SODM - Destination effects - Quarterly",
                          "SODM - Destination effects - Yearly","SODM - Destination and sector effects - Quarterly ")
names(source_labels) <- names_labels

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

# CReate the plot
p <- volatility_comparison_plot_tab %>% mutate(component = ifelse(component == 'Idio','Idiosyncratic', component)) %>%
                          ggplot() +
                          geom_line(aes(x = as.Date(period), y = as.numeric(value), group = component,
                                        linetype = component)) +  
                          geom_ribbon(aes(x = as.Date(period),
                                          ymin = as.numeric(value) - 1.96*as.numeric(std_err),
                                          ymax = as.numeric(value) + 1.96*as.numeric(std_err),
                                          group = component, color = component, fill = component,
                                          linetype = component), alpha = 0.5) +
                          theme_minimal() +
                          theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 18),
                                axis.text.y = element_text(size = 18),
                                panel.grid.major.y = element_line(color = 'azure2'),
                                panel.background = element_blank(),
                                axis.title.x = element_blank(),
                                axis.title.y = element_blank(),
                                legend.title = element_blank(),
                                legend.key= element_blank(),
                                legend.text = element_text(size = 22),
                                strip.text.x = element_text(
                                  size = 20, color = "black", face = "bold"
                                ),
                                strip.text.y = element_text(
                                  size = 20, color = "black", face = "bold"
                                )) +
                          scale_x_date(limits=c(min_d,max_d), date_labels = "%b %Y", breaks = scales::pretty_breaks(n = 10)) +
                          scale_color_grey(start = 0.4,
                                           end = 0.8) +
                          scale_fill_grey(start = 0.4,
                                          end = 0.8) +
                          facet_wrap(~source, ncol = 2, labeller = labeller(source = source_labels))
print(p)

# Run dev.off() to create the file!
dev.off()

##### Plot: Common growth rates components DFMs vs SODMs comparison boxplots (Figure 6) #####

str(matrix_list$full_dataset)

dates <-  seq(as.Date("1994-1-1"), as.Date("2017-10-1"), by = "quarter")

min_d <- as.Date("1994-1-1","%Y-%m-%d")
max_d <- as.Date("2017-10-31","%Y-%m-%d")

boxplot_data <- as.data.frame(rbind(cbind(source = 'DFM - Global and destination factors - Quarterly',
                                          component = 'Common',
                                          period = as.character(dates),
                                          colQuantiles(matrix_list$full_dataset$common, probs = c(0,.05,0.10,0.25,0.5,0.75,.9,.95,1), na.rm = T)),
                                    cbind(source = 'DFM - Global and destination factors - Quarterly',
                                          component = 'Idiosyncratic',
                                          period = as.character(dates),
                                          colQuantiles(matrix_list$full_dataset$idio, probs =c(0,.05,0.10,0.25,0.5,0.75,.9,.95,1), na.rm = T)),
                                    cbind(source = 'SODM - Destination and sector effects - Quarterly',
                                          component = 'Common',
                                          period = as.character(dates),
                                          colQuantiles(matrix_list$rep_destsec_quarterly$common, probs = c(0,.05,0.10,0.25,0.5,0.75,.9,.95,1), na.rm = T)),
                                    cbind(source = 'SODM - Destination and Sector effects - Quarterly',
                                          component = 'Idiosyncratic',
                                          period = as.character(dates),
                                          colQuantiles(matrix_list$rep_destsec_quarterly$idio, probs = c(0,.05,0.10,0.25,0.5,0.75,.9,.95,1), na.rm = T))))

colnames(boxplot_data)[4:12] <- c('MIN','Q5','Q10','Q25','Med','Q75','Q90','Q95','MAX')
rownames(boxplot_data) <- NULL

boxplot_data <- boxplot_data %>% mutate_at(c('MIN','Q5','Q10','Q25','Med','Q75','Q90','Q95','MAX'),as.numeric) %>%
                                 group_by(source, component, period) %>%
                                 mutate(IQR = Q25 - Q75) %>% 
                                 mutate(MIN_out = max(MIN, Q25 - 1.5*IQR)) %>%
                                 mutate(MAX_out = min(MAX, Q75 + 1.5*IQR)) 

# Call the pdf command to start the plot

pdf(file = here('results_analysis','figures','Fig_6_growth_rates_comparison.pdf'),   # The directory you want to save the file in
    width = 18, # The width of the plot in inches
    height = 8) # The height of the plot in inches

# Create the plot

p <- ggplot(boxplot_data) +
        geom_ribbon(aes(x = as.Date(period),
                        ymin = Q5, ymax = Q95,
                        group = component), alpha = 0.3) +
        geom_boxplot(aes(x = as.Date(period), 
                         ymin = Q10, lower = Q25, middle = Med, upper = Q75, ymax = Q90, fill = component,
                         group = interaction(period,component)), stat = 'identity',
                         position = position_dodge(width=0.9), color = 'black') +
        theme_minimal() +
        theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 18),
              axis.text.y = element_text(size = 18),
              panel.grid.major.y = element_line(color = 'azure2'),
              panel.grid.major.x = element_blank(),
              panel.background = element_blank(),
              axis.title.x = element_blank(),
              axis.title.y = element_blank(),
              legend.title = element_text(size = 22),
              legend.key= element_blank(),
              legend.text = element_text(size = 18),
              strip.text.x = element_text(size = 20, color = "black", face = "bold"),
              strip.text.y = element_text(size = 20, color = "black", face = "bold")) +
        scale_x_date(limits=c(min_d,max_d), date_labels = "%b %Y", breaks = scales::pretty_breaks(n = 10)) +
        scale_fill_grey(start = 0.5,
                        end = 1, name = 'Growth rates \ncomponent') +
        facet_wrap(component~source, ncol = 2,scales = 'free',
                   labeller = labeller(component = function(x) {rep("", length(x))}))
print(p)

# Run dev.off() to create the file!
dev.off()

##### Plot: Volatility components time series. Global, Destination and idiosyncratic components (Figure 8 and alternative versions) ######################

volatility_comparison_ts_tab <- vol_comparison_ts_tab %>% 
                                  filter((component %in% c('Global', 'Destination', 'Idio'))) %>% 
                                  filter(source == 'DFM Complete Dataset') 

# Call the pdf command to start the plot

pdf(file = here('results_analysis','figures','Fig_8_vol_components.pdf'),   # The directory you want to save the file in
    width = 18, # The width of the plot in inches
    height = 8) # The height of the plot in inches

# Create the plot
                                 
p <- volatility_comparison_ts_tab %>% mutate(component = ifelse(component == 'Idio','Idiosyncratic', component)) %>%
                                 mutate(component = factor(component, levels = c("Idiosyncratic","Global","Destination"))) %>%
                                  ggplot() +
                                 geom_ribbon(aes(x = as.Date(period),
                                                ymin = as.numeric(value) - 1.96*as.numeric(std_err),
                                                ymax = as.numeric(value) + 1.96*as.numeric(std_err),
                                                group = component, fill = component), color = 'black', alpha = 0.4) +
                                  geom_line(aes(x = as.Date(period), y = as.numeric(value), group = component,
                                                linetype = component), color = 'black') + 
                                  geom_point(aes(x = as.Date(period), y = as.numeric(value), group = component,
                                                 shape = component), cex = 3, color = 'black') +
                                  theme_minimal() +
                                  theme(plot.title = element_text(color="black", size=26, face="bold",
                                                                  hjust = 0.5),
                                        axis.text.x = element_text(angle = 45, hjust = 1, size = 18),
                                        axis.text.y = element_text(size = 18),
                                        panel.grid.major.y = element_line(color = 'azure2'),
                                        panel.background = element_blank(),
                                        axis.title.x = element_blank(),
                                        axis.title.y = element_blank(),
                                        legend.position = 'top',
                                        legend.title = element_blank(),
                                        legend.key= element_blank(),
                                        legend.text = element_text(size = 22),
                                        strip.text.x = element_text(
                                          size = 12, color = "black", face = "bold"
                                        ),
                                        strip.text.y = element_text(
                                          size = 12, color = "black", face = "bold"
                                        )) +
                                  ggtitle("Volatility components") +
                                  scale_x_date(limits=c(min_d,max_d), date_labels = "%b %Y",
                                               breaks = scales::pretty_breaks(n = 10)) +
                                  scale_fill_hue(c = 70, l = 60)

print(p) 
 
# Run dev.off() to create the file!
dev.off()

                                  
volatility_comparison_ts_tab %>% mutate(component = ifelse(component == 'Idio','Idiosyncratic', component)) %>%
                                  mutate(component = factor(component, levels = c("Idiosyncratic","Global","Destination"))) %>%
                                  ggplot() +
                                  geom_ribbon(aes(x = as.Date(period),
                                                  ymin = as.numeric(value) - 1.96*as.numeric(std_err),
                                                  ymax = as.numeric(value) + 1.96*as.numeric(std_err),
                                                  group = component, fill = component), color = 'black', alpha = 0.4) +
                                 geom_line(aes(x = as.Date(period), y = as.numeric(value), group = component,
                                                linetype = component), color = 'black') +
                                  theme_minimal() +
                                  theme(plot.title = element_text(color="black", size=26, face="bold",
                                                                  hjust = 0.5),
                                        axis.text.x = element_text(angle = 45, hjust = 1, size = 18),
                                        axis.text.y = element_text(size = 18),
                                        panel.grid.major.y = element_line(color = 'azure2'),
                                        panel.background = element_blank(),
                                        axis.title.x = element_blank(),
                                        axis.title.y = element_blank(),
                                        legend.position = 'top',
                                        legend.title = element_blank(),
                                        legend.key= element_blank(),
                                        legend.text = element_text(size = 22),
                                        strip.text.x = element_text(
                                          size = 12, color = "black", face = "bold"
                                        ),
                                        strip.text.y = element_text(
                                          size = 12, color = "black", face = "bold"
                                        ))  +
                                  ggtitle("Volatility components") +
                                  scale_x_date(limits=c(min_d,max_d), date_labels = "%b %Y",
                                               breaks = scales::pretty_breaks(n = 10)) +
                                  scale_fill_grey(start = 0, end = .9)
                                

############### Summary table with point estimates of the volatility over the time span (Table 1) ######################

agg_vol_from_samples <- agg_vol_tab_samples %>% gather(agg_type, VALUE, 2:ncol(.)) %>% group_by(agg_type) %>%
                         mutate(AVG = mean(VALUE))  %>%
                         mutate(bts_cmax90 = bootstrap_CI(VALUE,0.95)) %>%
                         mutate(bts_cmin90 = bootstrap_CI(VALUE,0.05)) %>%
                         mutate(bts_cmax99 = bootstrap_CI(VALUE,0.99)) %>%
                         mutate(bts_cmin99 = bootstrap_CI(VALUE,0.01)) %>%
                         mutate(MIN =  min(VALUE)) %>%
                         mutate(MAX = max(VALUE)) %>%
                         mutate(CI=ci_fun(VALUE))%>%
                         mutate(p_CI = CI/AVG) %>%
                         ungroup() 

agg_vol_from_samples <- agg_vol_from_samples[,-c(1,3)] %>% unique()

data_selection_table <- agg_vol_from_samples %>% mutate(agg_tech = substr(agg_type,1,2)) %>% 
                         mutate(comp = substr(agg_type,8,nchar(agg_type))) %>%
                         filter(comp %in% c('common_na','dest_na','glob_na','idio')) %>%  
                         mutate(across(where(is.numeric), function(x) round(x, 4))) 
setDT(data_selection_table)
output_table <- data.frame()

components <- c('common_na','glob_na','dest_na','idio')
pretty_names <- c('Common','Global','Destination','Idiosyncratic')
for(cmp in components)
{  
index <- which(components == cmp)  
## Fill the common component results
dynamic_w_row <- data_selection_table[comp == cmp & agg_tech == 'M1']
constant_w_row <- data_selection_table[comp == cmp & agg_tech == 'M3']
print_conf_interval <- function(x){paste0('(',format(x$AVG + x$bts_cmin99,nsmall = 4),',',
                                          format(x$AVG + x$bts_cmax99,nsmall = 4),')')}
output_table['Components',2*index] <- output_table['Components',2*index-1] <- pretty_names[index]
output_table['Aggregation method',2*index-1] <- 'Dynamic weights'
output_table['Aggregation method',2*index] <- 'Constant weights'
output_table['DFM (sampled)',2*index-1] <- dynamic_w_row$AVG
output_table['DFM (sampled)',2*index] <- constant_w_row$AVG
output_table['CI (sampled)',2*index-1] <- print_conf_interval(dynamic_w_row)
output_table['CI (sampled)',2*index] <- print_conf_interval(constant_w_row)
}  

### Filling with the reference values for SODM
repl_data_tab <- replication_vol_tab %>% mutate(across(2:ncol(.),as.numeric)) %>%
                                         mutate(across(where(is.numeric), function(x) round(x, 4))) %>%
                                         mutate(across(where(is.numeric), function(x) format(x, nsmall = 4))) 

for(model in repl_data_tab$SIM)
{  
output_table[model,1] <- repl_data_tab[SIM == model]$M1_agg_common
output_table[model,2] <- repl_data_tab[SIM == model]$M3_agg_common
output_table[model,7] <- repl_data_tab[SIM == model]$M1_agg_idio
output_table[model,8] <- repl_data_tab[SIM == model]$M3_agg_idio
}


print(xtable(t(output_table), type = "latex"))


fileConn <- file(here('results_analysis','tables',"Tab_1_models_comparison.txt"))
writeLines(print(xtable(t(output_table)), type = "latex"), fileConn)
close(fileConn)


#########################################################
########### Plots and tables for Section 4.2 ############
#########################################################

#####  Calculate pvalues for block specific volatilities  #####
cnt_vols_long <- cnt_vols_tab %>% gather(TYPE,VOL,4:ncol(cnt_vols_tab)) %>% 
                                  mutate(TYPE = substr(TYPE,5,nchar(TYPE))) %>%
                                  mutate(VOL = sqrt(as.numeric(VOL)))

cnt_vols_long$VOL <- as.numeric(cnt_vols_long$VOL)
cnt_vols_long$FLOWS <- as.numeric(cnt_vols_long$FLOWS)


cnt_vols_pval <- cnt_vols_long %>% filter(SIM != 0) %>% 
                                  group_by(PYOD,TYPE) %>%  mutate(AVG_FLOWS = mean(FLOWS)) %>%
                                  mutate(AVG_VOL = mean(VOL))  %>%
                                  mutate(cmax90 = bootstrap_CI(VOL,0.95)) %>%
                                  mutate(cmin90 = bootstrap_CI(VOL,0.05)) %>%
                                  mutate(cmax99 = bootstrap_CI(VOL,0.99)) %>%
                                  mutate(cmin99 = bootstrap_CI(VOL,0.01)) %>%
                                  mutate(MIN =  min(VOL)) %>%
                                  mutate(MAX = max(VOL)) %>%
                                  mutate(CI=ci_fun(VOL))%>%
                                  mutate(p_CI = CI/AVG_VOL) %>%
                                  ungroup() 

setDT(cnt_vols_pval)

cnt_vols_pval <- unique(cnt_vols_pval[,-c('SIM','FLOWS','VOL'), with = F]) 

###### Prepare the dataset for gravity-like analysys for the volatilities associated to each destination ######

load(here::here('results_analysis','geo_data.Rdata')) ### import geo data downloaded from 
setDT(geo_data)

setkey(geo_data,PYOD)
setkey(cnt_vols_pval, PYOD)

cnt_vols_geo <- geo_data[cnt_vols_pval]
cnt_vols_geo$eu_d[is.na(cnt_vols_geo$eu_d)] <- 0


cnt_vols_geo$colony[cnt_vols_geo$PYOD == 'US'] <- 0
cnt_vols_geo <- within(cnt_vols_geo, fta_bb[col45 == 1] <- 6) 
cnt_vols_geo <- within(cnt_vols_geo, fta_bb[(col45 != 1)&(colony == 1)] <- 5)
cnt_vols_geo$eu_col <- cnt_vols_geo$eu_d
cnt_vols_geo <- within(cnt_vols_geo, eu_col[col45 == 1] <- 2) 
cnt_vols_geo <- within(cnt_vols_geo, eu_col[(col45 != 1)&(colony == 1)] <- 3)

setDT(cnt_vols_geo)
cnt_vols_geo$fta_bb[is.na(cnt_vols_geo$fta_bb)] <- 0
cnt_vols_geo$fta_bb <- factor(cnt_vols_geo$fta_bb)
cnt_vols_geo$eu_col <- factor(cnt_vols_geo$eu_col)
cnt_vols_geo$eu_d <- factor(cnt_vols_geo$eu_d)

pacman::p_load(plyr)
cnt_vols_geo <- cnt_vols_geo %>% mutate(fta_verbose = revalue(fta_bb, c('0'='None','1'='FTA', '2' = 'Custom union', '3' = 'Common market', '4' = 'Common economy', '5' = 'Ex-colony', '6'='Ex-colony (after 1945)'))) %>% 
  mutate(eu_col_verbose = revalue(eu_col, c('0'='Extra-EU','1'='Intra-EU', '2' = 'Ex colony (after 1945)', '3' = 'Ex colony'))) %>%
  mutate(eu_verbose = revalue(eu_d, c('0'='Extra-EU','1'='Intra-EU')))
detach("package:plyr", unload=TRUE)    


############### Plot with destination-specific volatilities vs GDP (Figure 9) ######################

pdf(file = here('results_analysis','figures','Fig_9_GDP_vs_vol.pdf'),   # The directory you want to save the file in
    width = 18, # The width of the plot in inches
    height = 8) # The height of the plot in inches

# Create the plot
p <- cnt_vols_geo %>% filter(TYPE == 'CONST_WGT_RES') %>% 
                  ggplot() +
                  geom_point(aes(x = gdp_d/10^9, y = AVG_VOL, size = 100, color = eu_col_verbose, shape = eu_col_verbose), alpha = 0.2, cex = 14)  +
                  geom_text(aes(x=gdp_d/10^9,y=AVG_VOL, color=eu_col_verbose, label = PYOD, size = 14))+
                  scale_x_continuous(trans ='log', labels = fmt_decimals2(2)) +  
                  scale_y_continuous(trans = 'log', labels = fmt_decimals2(2)) +
                  scale_color_hue(l = 40, c=50) + 
                  scale_shape_manual(values = c(15,17,18,19)) +
                  theme_minimal()+ 
                  theme() +
                  theme(axis.text.x = element_text(angle = 45, hjust = 1, 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 = '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")) +
                  guides(color = guide_legend(nrow = 1)) +
                  scale_size_continuous(guide = F)+
                  labs(y = 'Destination-specific volatility', x = 'Destination GDP (Billion $)', title = '') +
                  facet_wrap(~eu_verbose)
print(p)

# Run dev.off() to create the file!
dev.off()

#### Plot with destination-specific volatilities vs GDP (Robustness check with time-varying weights aggregation) ####

cnt_vols_geo %>% filter(TYPE == 'DYN_WGT_RES') %>% 
                  ggplot() +
                  geom_point(aes(x = gdp_d/10^9, y = AVG_VOL, size = 100, color = eu_col_verbose, shape = eu_col_verbose), alpha = 0.2, cex = 14)  +
                  geom_text(aes(x=gdp_d/10^9,y=AVG_VOL, color=eu_col_verbose, label = PYOD, size = 14))+
                  scale_x_continuous(trans ='log', labels = fmt_decimals2(2)) +  
                  scale_y_continuous(trans = 'log', labels = fmt_decimals2(2)) +
                  scale_color_hue(l = 40, c=50) + 
                  scale_shape_manual(values = c(15,17,18,19)) +
                  theme_minimal()+ 
                  theme() +
                  theme(axis.text.x = element_text(angle = 45, hjust = 1, 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 = '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")) +
                  guides(color = guide_legend(nrow = 1)) +
                  scale_size_continuous(guide = F)+
                  labs(y = 'Destination-specific volatility', x = 'Destination GDP (Billion $)', title = '') +
                  facet_wrap(~eu_verbose)


#### Gravity equation regression (Table 2) ####

cnt_vols_geo_reg <- cnt_vols_geo %>% group_by(TYPE) %>% mutate(q_flows = cut(log(AVG_FLOWS),quantile(log(AVG_FLOWS), probs = 1:5/5))) %>% ungroup()
unique(cnt_vols_geo_reg[,c('PYOD','q_flows')])

ty <- unique(cnt_vols_geo_reg$TYPE)[6]

d1 <- lm(log(AVG_VOL) ~ log(distw) + log(gdp_d), data=filter(cnt_vols_geo_reg, TYPE== ty))
d2 <- lm(log(AVG_VOL) ~ log(distw) + log(gdp_d) + factor(eu_col_verbose), data=filter(cnt_vols_geo_reg, TYPE== ty))
d3 <- lm(log(AVG_VOL) ~ log(distw) + log(gdp_d) + factor(fta_verbose), data=filter(cnt_vols_geo_reg, TYPE== ty))

#tab_model(d1,d2,d3)


fileConn <- file(here('results_analysis','tables',"Tab_2_gravity_regression.txt"))
writeLines(print(stargazer::stargazer(d1,d2,d3)), fileConn)
close(fileConn)


##### Saving main results in R files ################

save(cnt_vols_tab, file = here::here('results_analysis','partial_outputs',"cnt_vols_tab.RData"))
save(fac_list, file = here::here('results_analysis','partial_outputs','fac_list.RData'))
save(agg_cor_list, file = here::here('results_analysis','partial_outputs',"agg_cor_list.RData" ))
save(agg_vol_tab_samples, file= here::here('results_analysis','partial_outputs',"agg_vol_tab_samples.RData"))
save(agg_vol_tab_complete_dataset, file= here::here('results_analysis','partial_outputs',"agg_vol_tab_complete_dataset.RData"))

