#### Auxiliary packages #####

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

pacman::p_load(here)
pacman::p_load(data.table)
pacman::p_load(tidyr)
pacman::p_load(haven)

f_nan = function(DT) 
{
  for (i in names(DT))
    DT[get(i)==0, (i):=NA]
}

### impose country selection
load(here('data_raw','countries.RData'))

###

data <- fread(here::here('data_raw','export_val_country_quart_dataset_1993-2017.csv'))
#data <- data[,-1]

export_cols <- colnames(data)[3:ncol(data)]
data[,(export_cols) := lapply(.SD, function(x) as.numeric(x)), .SDcols = export_cols]


## Import also the sectors for the integration with the dataset
sector_siren <- fread(here::here('data_raw','sectors_siren.csv'))

f_nan(data)

mat_base <- as.matrix(data[,3:ncol(data)]) ## extract the numeric matrix
gr_mat <- log(1+mat_base[,5:ncol(mat_base)]) - log(1+mat_base[,1:(ncol(mat_base)-4)]) ## calculate GR
gr_base <- setDT(as.data.frame(gr_mat))  ## convert to data.table

data <- data[rowSums(!is.na(gr_base), na.rm = T)>1] ## remove empty obs in GR from original dataset
gr_base <- gr_base[rowSums(!is.na(gr_base), na.rm = T)>1] ## remove empty obs in GR from GR dataset
gr_base <- cbind(data[,1:2, with = F],gr_base) ## insert siren pyod indicators


q_gr_datasets <- list()
periods <- colnames(gr_base)[3:ncol(gr_base)]
for(y in 1:length(periods))
{
  if(y %% 4 == 1) 
  {periods_vec <- periods[((0:3)+y)]  
  aux_mat <- gr_base[,c('SIREN','PYOD',periods), with = F]  
  long_mat <- tidyr::gather(aux_mat, period, gr_capays, 3:ncol(aux_mat))
  setDT(long_mat)
  q_gr_datasets[[y]] <- long_mat[!is.na(gr_capays)]}
}  
q_gr_dataset <- do.call(rbind, q_gr_datasets)
q_gr_dataset[,c('year','quarter'):= list(as.numeric(substring(period,1,4)),as.numeric(substring(period,6,7)))]
q_gr_dataset[, period_numeric := year + (quarter-1)*0.25]

q_capays_datasets <- list()
periods <- colnames(data)[3:ncol(data)]
for(y in 1:length(periods))
{
  if(y %% 4 == 1) 
  {periods_vec <- periods[((0:3)+y)]  
  aux_mat <- data[,c('SIREN','PYOD',periods_vec), with = F]  
  long_mat <- tidyr::gather(aux_mat, period, capays, 3:ncol(aux_mat))
  setDT(long_mat)
  q_capays_datasets[[y]] <- long_mat[!is.na(capays)]}
}  
q_capays_dataset <- do.call(rbind, q_capays_datasets)
q_capays_dataset[,c('year','quarter'):= list(as.numeric(substring(period,1,4)),
                                         as.numeric(substring(period,6,7)))]
q_capays_dataset[, period_numeric := year + (quarter-1)*0.25]


setkey(q_gr_dataset,SIREN,PYOD,year,quarter, period, period_numeric)
setkey(q_capays_dataset,SIREN,PYOD,year,quarter, period, period_numeric)


long_panel_quarterly <- q_gr_dataset[q_capays_dataset]

long_panel_quarterly[,n := .N, by = .(SIREN,PYOD)]
long_panel_quarterly <- long_panel_quarterly[n>4]
long_panel_quarterly[,lag_capays := c(NA,NA,NA,NA,capays[-c(.N,.N-1,.N-2, .N-3)]), keyby = .(SIREN,PYOD)]

columns <- colnames(long_panel_quarterly)
colnames(long_panel_quarterly)[which(columns == 'SIREN')] <- 'siren'
colnames(long_panel_quarterly)[which(columns == 'PYOD')] <- 'pays'

colnames(sector_siren)[which(columns == 'SIREN')] <- 'siren'
setkey(sector_siren,siren)
setkey(long_panel_quarterly, siren)

long_panel_quarterly <- sector_siren[long_panel_quarterly]

long_panel_quarterly <- long_panel_quarterly[pays %in% countries]

long_panel_quarterly[,out_gr := ifelse(capays/lag_capays > 2 | capays/lag_capays < 0.5,1,0)]
colnames(long_panel_quarterly)[sapply(colnames(long_panel_quarterly),function(x) grepl('sec',x, fixed = TRUE))] <- 'main_sector'
colnames(long_panel_quarterly)[sapply(colnames(long_panel_quarterly),function(x) grepl('NC2',x, fixed = TRUE))] <- 'main_sector'


#fwrite(long_panel_quarterly, here::here('DiGiovanni_et_al_replication','Input','long_panel_quarterly.csv'))
write_dta(long_panel_quarterly, here::here('DiGiovanni_et_al_replication','Input','long_panel_quarterly.dta'))

######## Work on yearly export flows #########

dataset_new <- data[,1:2]
colnames(dataset_new)[1:2] <- colnames(data)[1:2]
n_years <- floor((ncol(data)-2)/4)
for(s in 1:n_years)
{
  dataset_new <- cbind(dataset_new,rowSums(data[,2+(4*(s-1)+1):(4*(s-1)+4),with = F], na.rm = T))
  colnames(dataset_new)[2+s] <- substr(colnames(data)[2+(4*(s-1)+1)],1,4)
} 

f_nan(dataset_new) ## substituting NAs with zeros (at this stage they are equivalent)

mat_base <- as.matrix(dataset_new[,3:ncol(dataset_new)]) ## extract the numeric matrix
gr_mat <- log(1+mat_base[,2:ncol(mat_base)]) - log(1+mat_base[,1:(ncol(mat_base)-1)]) ## calculate GR
gr_base <- setDT(as.data.frame(gr_mat))  ## convert to data.table

dataset_new <- dataset_new[rowSums(!is.na(gr_base), na.rm = T)>1] ## remove empty obs in GR from original dataset
gr_base <- gr_base[rowSums(!is.na(gr_base), na.rm = T)>1] ## remove empty obs in GR from GR dataset
gr_base <- cbind(dataset_new[,1:2, with = F],gr_base) ## insert siren pyod indicators

gr_datasets <- list()
years <- colnames(gr_base)[3:ncol(gr_base)]
for(y in years)
{
y <- as.numeric(y)
if(as.numeric(y) %% 2 == 0) 
{years_vec <- as.character((0:1)+y)  
aux_mat <- gr_base[,c('SIREN','PYOD',years_vec), with = F]  
long_mat <- tidyr::gather(aux_mat, year, gr_capays, 3:ncol(aux_mat))
setDT(long_mat)
gr_datasets[[y]] <- long_mat[!is.na(gr_capays)]}
}  
gr_dataset <- do.call(rbind, gr_datasets)


dataset_pays_list <- NULL
years <- colnames(dataset_new)[3:ncol(dataset_new)]
for(y in years)
{
  y <- as.numeric(y)
  if(as.numeric(y) %% 5 == 3) 
  {years_vec <- as.character((0:4)+y)  
  aux_mat <- dataset_new[,c('SIREN','PYOD',years_vec), with = F]  
  long_mat <- tidyr::gather(aux_mat, year, capays, 3:ncol(aux_mat))
  setDT(long_mat)
  dataset_pays_list[[y]] <- long_mat[!is.na(capays)]}
} 
dataset_pays <- do.call(rbind, dataset_pays_list)

setkey(dataset_pays,SIREN,PYOD,year)
setkey(gr_dataset,SIREN,PYOD,year)

long_panel_yearly <- gr_dataset[dataset_pays]


long_panel_yearly[,n := .N, by = .(SIREN,PYOD)]
long_panel_yearly <- long_panel_yearly[n>2]
long_panel_yearly[,lag_capays := c(NA,capays[-.N]), keyby = .(SIREN,PYOD)]

columns <- colnames(long_panel_yearly)
colnames(long_panel_yearly)[which(columns == 'SIREN')] <- 'siren'
colnames(long_panel_yearly)[which(columns == 'PYOD')] <- 'pays'

colnames(sector_siren)[which(columns == 'SIREN')] <- 'siren'
setkey(sector_siren,siren)
setkey(long_panel_yearly, siren)

long_panel_yearly <- sector_siren[long_panel_yearly]

long_panel_yearly <- long_panel_yearly[pays %in% countries]

long_panel_yearly[,out_gr := ifelse(capays/lag_capays > 2 | capays/lag_capays < 0.5,1,0)]
colnames(long_panel_yearly)[sapply(colnames(long_panel_yearly),function(x) grepl('sec',x, fixed = TRUE))] <- 'main_sector'
colnames(long_panel_yearly)[sapply(colnames(long_panel_yearly),function(x) grepl('NC2',x, fixed = TRUE))] <- 'main_sector'


#fwrite(long_panel_yearly, here::here('DiGiovanni_et_al_replication','Input','long_panel_yearly.csv'))
write_dta(long_panel_yearly, here::here('DiGiovanni_et_al_replication','Input','long_panel_yearly.dta'))
