################################################################################
###### Main script for the preliminary analysis of the dataset properties ######
################################################################################

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

#### Auxiliary packages ####
pacman::p_load(here)
pacman::p_load(data.table)
pacman::p_load(tidyr)
pacman::p_load(dplyr)
pacman::p_load(ggplot2)
pacman::p_load(matrixStasts)

#### Import the dataset from the data_raw folder ####
dataframe <- fread(here::here('data_raw','export_val_country_quart_dataset_1993-2017.csv'), )
#dataframe[,V1 := NULL] ### to remove the varibale associated to row numbers if uncomment if row numbers are imported

## Define a function converting the NAs to zeros
f_zero = function(DT) 
{
  for (i in names(DT))
    DT[is.na(get(i)), (i):=0]
}
## Define a function converting the zeros to NAs
f_nan = function(DT) 
{
  for (i in names(DT))
    DT[get(i)==0, (i):=NA]
}

# Set all missings od principal dataset to zeros
f_zero(dataframe)

export_value_cols <- colnames(dataframe)[4:ncol(dataframe)]


##### Counting NAs in the dataset #####
f_nan(dataframe)

## Calculating the NA share in the dataset and printing the outcome
na_share <- sum(rowSums(is.na(dataframe[,3:ncol(dataframe)])))/(100*nrow(dataframe))
print(paste0('Percentage of missing values in the orginal dataset: ',round(na_share*100,2),'%'))

## Computing the NA shares per quarter and plotting the time series

col_shares <- colSums(is.na(dataframe[,3:ncol(dataframe)]))/(nrow(dataframe))
na_plot_data <- as.data.frame(cbind(Date = names(col_shares), Value = col_shares))

## Fromatting the Date type for plotting
pacman::p_load(stringr)
na_plot_data$Date <- str_replace_all(na_plot_data$Date,".01$","-01-01")
na_plot_data$Date <- str_replace_all(na_plot_data$Date,".04$","-10-01")
na_plot_data$Date <- str_replace_all(na_plot_data$Date,".02$","-04-01")
na_plot_data$Date <- str_replace_all(na_plot_data$Date,".03$","-07-01")

na_plot_data$Date <- as.Date(na_plot_data$Date,origin = "1970-01-01")

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

## Producing the plot of the time series of the NA shares per quarter

pdf(file = here('data_preparation_sampling','data_summary','NA_vs_time_plot.pdf'),   # The directory you want to save the file in
    width = 14, # The width of the plot in inches
    height = 8) # The height of the plot in inches

# Create the plot
p <- na_plot_data %>%
            ggplot() +
            geom_line(aes(x = Date, y = as.numeric(Value), color = 'red'), size = 1.1) +
            theme_light() +
            theme(axis.text.x = element_text(angle = 45, hjust = 1),
                  panel.grid.minor.y = element_blank(),
                  legend.position = 'none')+
            scale_color_hue(l=40, c=35) +
            labs(y = '', x = '', title = '') +
            scale_x_date(limits = c(min_d,max_d),date_labels = "%b/%Y", breaks = scales::pretty_breaks(n = 14)) 
print(p)
dev.off()

## Computing the empirical distributions of number of flows with
vec_na_num <- rowSums(is.na(dataframe[,3:ncol(dataframe)])) ### computing the number of NA per flow
vec_na_num <- as.data.frame(cbind(x = 1:length(vec_na_num), N_NA  = vec_na_num))
setDT(vec_na_num)
n_flow_dist <- vec_na_num[,.N, keyby = c('N_NA')]
colnames(n_flow_dist)[2] <- 'N_flows'
tot_flows <- nrow(vec_na_num)
rm(vec_na_num)


## Plotting the distribution on a log-log scale
n_flow_dist[N_NA <100] %>%
                      ggplot() +
                      geom_point(aes(x = 1-(N_NA/100), y = N_flows*100/tot_flows))  +
                      scale_y_continuous(trans = 'log') +
                      scale_x_continuous(trans = 'log') +
                      scale_color_hue(l=40, c=35) +
                      theme_minimal() +
                      labs(y = '', x = '', title = '')


############ Computing few summary statistics of the dataset ##############

source(here::here('data_preparation_sampling','data_summary_function.R'))

data_summary_function(dataframe, 'data_summary')
tab_sum <- fread(here::here('data_preparation_sampling','data_summary','data_summary_by_siren.csv'))[,-1]
statistics_and_volatility(tab_sum, 'data_summary')


############ Destination selection and sampling ################

source(here::here('data_preparation_sampling','destinations_selection_and_sampling_function.R'))

# Define the list of parameteres for the country selection procedure and the sampling
filt_par <- list(flow_filter = 100,        # Minimum number of series for destinations to be included (flow_filter = 1000)
                 share_of_exp_value_filter = 0.012,   # Minimum share of export in a given quarter to be included (share_of_exp_value_filter = 0.012)
                 num_firms = 5000,  # Firms selected in each sample (num_firms = 100000)
                 num_samples = 20)   # Number of samples

sampling  <- T
base_sel <- destinations_selection_and_sampling(dataframe,filt_par, sampling)

