# КОД ДЛЯ ВОСПРОИЗВОДСТВА РЕЗУЛЬТАТОВ СТАТЬИ:
# Гнидченко А.А. Концентрация мировой торговли и сегрегация товарных 
# рынков // Журнал Новой экономической ассоциации, 2025, №1 (66).

# <используемая кодировка: UTF-8> #


# 1) Подготовительный этап ####
## 1.1) Автоматическая установка директории ####
bx <- "box"
bx_new <- bx[!(bx %in% installed.packages()[,"Package"])]
if(length(bx_new)) install.packages(bx_new)
setwd(box::file())
getwd() # проверка


## 1.2) Загрузка пакетов ####
packs <- c(
  "data.table","dplyr","Hmisc","grid","gridExtra","ggplot2","MASS","car",
  "Cairo","factoextra","gtools","openxlsx","glue","stats","reshape2")
newpacks <- packs[!(packs %in% installed.packages()[,"Package"])]
if(length(newpacks)) install.packages(newpacks)
lapply(packs, library, character.only=T)[[1]]


## 1.3) Прописывание пользовательских функций ####
`%!in%` = Negate(`%in%`)
zv <- function(x){ifelse(x<0.01,"***",ifelse(x<0.05,"**",ifelse(x<0.1,"*","")))}
NAP <- function(x) { x[is.na(x)] <- ""; x }
HHI <- function(data) {
  data <- data.frame(as.numeric(as.character(data)))
  res <- sum((data/sum(data, na.rm=T))^2, na.rm=T); res
}
sumNA <- function(value,threshold) {
  value=as.numeric(value)
  if(missing(threshold)) {
    sum(value,na.rm = TRUE)
  } else {
    ifelse(threshold>1 | threshold<0,
           stop("threshold should vary between 0 and 1"),
           ifelse(threshold>=sum(is.na(value)==TRUE)/length(value),
                  sum(value,na.rm = TRUE),sum(value,na.rm = FALSE)))
  }
}
medianNA <- function(value,threshold) {
  if(missing(threshold)) {
    median(value,na.rm = TRUE)
  } else {
    ifelse(threshold>1 | threshold<0,
           stop("threshold should vary between 0 and 1"),
           ifelse(threshold>=sum(is.na(value)==TRUE)/length(value),
                  median(value,na.rm = TRUE),median(value,na.rm = FALSE)))
  }
}
meanNA <- function(value,threshold) {
  if(missing(threshold)) {
    mean(value,na.rm = TRUE)
  } else {
    ifelse(threshold>1 | threshold<0,
           stop("threshold should vary between 0 and 1"),
           ifelse(threshold>=sum(is.na(value)==TRUE)/length(value),
                  mean(value,na.rm = TRUE),mean(value,na.rm = FALSE)))
  }
}
w.perc <- function (x, weights, probs, normwt = F, na.rm,
                    type = c("quantile", "(i-1)/(n-1)", "i/(n+1)", "i/n")) 
{
  if (!length(weights)) 
    return(quantile(x, probs = probs, na.rm = na.rm))
  type <- match.arg(type)
  if (any(probs < 0 | probs > 1)) 
    stop("Probabilities must be between 0 and 1 inclusive")
  nams <- paste(format(round(probs * 100, if (
    length(probs) > 1) 2 - log10(diff(range(probs))) else 2)), "%", sep = "")
  if (type == "quantile") {
    library("Hmisc")
    w <- wtd.table(x, weights, na.rm = na.rm, normwt = normwt, 
                   type = "list")
    x <- w$x
    wts <- w$sum.of.weights
    n <- sum(wts)
    order <- 1 + (n - 1) * probs
    low <- pmax(floor(order), 1)
    high <- pmin(low + 1, n)
    order <- order%%1
    allq <- tryCatch(
      {
        approx(cumsum(wts), x, xout = c(low, high), method = "constant", 
               f = 1, rule = 2)$y 
      },
      error=function(cond) {
        message("Some minor errors executed:")
        message(cond)
        message("")
        return(NA)
      },
      finally={
      }
    )
    k <- length(probs)
    quantiles <- (1 - order) * allq[1:k] + order * allq[-(1:k)]
    quantiles <- if (length(quantiles)==k) {
      quantiles
    } else {
      vector(length=k)*NA
    }
    names(quantiles) <- nams
    return(quantiles)
  }
  w <- wtd.Ecdf(x, weights, na.rm = na.rm, type = type, normwt = normwt)
  structure(approx(w$ecdf, w$x, xout = probs, rule = 2)$y, 
            names = nams)
}
theme_GOST_base <- function (base_size = 14, base_family = "serif") {
  theme_bw(base_size = base_size, 
           base_family = base_family) +
    theme(
      plot.title = element_text(hjust = 0.5, size = 16,
                                margin = margin(t = 0, b = 10)),
      axis.text = element_text(color = "black"),
      plot.margin = margin(0.5, 0.2, 0.25, 0.25, "cm")
    )
}



# 2) Подготовка данных ####
## 2.1) Загрузка данных с сайта CEPII ####
getOption('timeout')
options(timeout=100)
temp <- tempfile()
download.file("http://www.cepii.fr/DATA_DOWNLOAD/baci/data/BACI_HS17_V202401.zip",temp)
CepiiData <- read.csv(unz(temp, "BACI_HS17_Y2022_V202401.csv"))
do.call(file.remove, list(temp))
unlink(temp)
gc()

year <- CepiiData$t[1]


## 2.2) Первичная обработка данных ####
CepiiData <- CepiiData[, -which(names(CepiiData) %in% c("t","q"))]
CepiiData$k[nchar(CepiiData$k)==5] <- paste0("0",CepiiData$k[nchar(CepiiData$k)==5])
CepiiData$i <- as.character(CepiiData$i)
CepiiData$j <- as.character(CepiiData$j)
CepiiData$c <- substr(CepiiData$k,1,4)
gc()


## 2.3) Формирование данных на уровне 4 знаков ТН ВЭД ####
Cepii4 <- data.table::dcast(as.data.table(CepiiData), i+j+c~"v", sum, value.var=c("v"))
names(Cepii4)[ncol(Cepii4)] <- "v"
rm(CepiiData)
gc()



# 3) Расчеты ####
## 3.1) Расчет мирового объема торговли по товару ####
TOT <- Cepii4 %>% group_by(c) %>% summarise(v=sum(v)) %>% as.data.frame()
TOT$v <- TOT$v/1000000 # млрд долл.


## 3.2) Расчет базовых показателей ####
# (этот блок считается несколько минут)
setDT(Cepii4)[, V_m := sumNA(v), by = .(j, c)] # объем импорта по товару c из всех стран
setDT(Cepii4)[, V_x := sumNA(v), by = .(i, c)] # объем экспорта по товару c во все страны
setDT(Cepii4)[, hi_m := HHI(v), by = .(j, c)] # HHI для импортера товара c по странам-экспортерам
setDT(Cepii4)[, hi_x := HHI(v), by = .(i, c)] # HHI для экспортера товара c по странам-импортерам
gc()


## 3.3) Рис. 1 ####
options(scipen=7)
tov_exmpl <- c("8601","2606")
tit_exmpl <- c("Локомотивы железнодорожные (ТН ВЭД 8601)",
               "Руды и концентраты алюминиевые (ТН ВЭД 2606)")
leg_pos <- c("none","bottom")
list_pic <- list()
i=1
while(i<length(tov_exmpl)+1){
  tov <- tov_exmpl[i]
  test <- subset(Cepii4, c == tov)
  test2 <- gtools::mixedsort(unique(c(test$i, test$j)))
  test2 <- as.data.frame(test2)
  names(test2) <- "country"
  test2$hi_m <- test$hi_m[match(test2$country, test$j)]
  test2$hi_x <- test$hi_x[match(test2$country, test$i)]
  test2$V_m <- test$V_m[match(test2$country, test$j)]
  test2$V_x <- test$V_x[match(test2$country, test$i)]
  medM <- medianNA(test2$hi_m) # медиана | импорт
  medX <- medianNA(test2$hi_x) # медиана | экспорт
  meanM <- meanNA(test2$hi_m) # средняя | импорт
  meanX <- meanNA(test2$hi_x) # средняя | экспорт
  wmedM <- w.perc(test2$hi_m, test2$V_m, 0.5, na.rm=T) # взвешенная медиана | импорт
  wmedX <- w.perc(test2$hi_x, test2$V_x, 0.5, na.rm=T) # взвешенная медиана | экспорт
  wmeanM <- sumNA(test2$hi_m*test2$V_m)/sumNA(test2$V_m)
  wmeanX <- sumNA(test2$hi_x*test2$V_x)/sumNA(test2$V_x)
  nM <- length(test2$hi_m[is.na(test2$hi_m)==F])
  nX <- length(test2$hi_x[is.na(test2$hi_x)==F])
  wnM <- 1/HHI(test2$V_m) # взвешенное количество импортеров (1/HHI)
  wnX <- 1/HHI(test2$V_x) # взвешенное количество экспортеров (1/HHI)
  
  test3 <- reshape2::melt(test2)
  test3$variable <- as.character(test3$variable)
  test3$flow <- substr(test3$variable,(nchar(test3$variable)),nchar(test3$variable))
  test3$var <- gsub("_.","",test3$variable)
  test3 <- reshape2::dcast(test3, country+flow~var, sum, value.var="value")
  test3$flow[test3$flow=="m"] <- "концентрация поставщиков"
  test3$flow[test3$flow=="x"] <- "концентрация покупателей"
  
  pic_test <- ggplot(test3) + geom_density(aes(hi, weight=V, col=flow, linewidth=flow)) + 
    scale_color_manual(values=c("grey50","black")) + scale_linewidth_manual(values=c(1,0.75)) + 
    geom_vline(xintercept=wmeanM, col="grey50", linetype="dashed") + 
    geom_vline(xintercept=wmeanX, col="black", linetype="dashed") +
    xlab("Индекс Херфиндаля-Хиршмена") + ylab("Плотность распределения \n(взвешенная торговлей)") + 
    theme_GOST_base() + ggtitle(tit_exmpl[i]) + scale_x_continuous(limits=c(0,1), breaks=seq(0,1,0.2)) +
    theme(legend.position=leg_pos[i], legend.title=element_blank(), legend.text=element_text(size=13))
  list_pic[[i]] <- pic_test
  i=i+1
}
png('./pic1.png', 
    width = 6, height = 8, units = 'in', res = 500)
gridExtra::grid.arrange(list_pic[[1]], list_pic[[2]], heights=c(3.7,4.3))
dev.off()
CairoPDF("./pic1.pdf", width = 6.5, height = 8.5)
gridExtra::grid.arrange(list_pic[[1]], list_pic[[2]], heights=c(3.7,4.3))
dev.off()
###


## 3.4) Расчетный цикл по всем товарам ####
TovList <- gtools::mixedsort(unique(Cepii4$c))
DF <- data.frame(matrix(ncol=13,nrow=length(TovList)))
names(DF) <- c("TNVED","medM","medX","wmedM","wmedX","meanM","meanX","wmeanM","wmeanX","nM","nX","wnM","wnX")
i=1
while(i<length(TovList)+1){
  Tov <- TovList[i]
  test <- subset(Cepii4, c==Tov)
  test2 <- gtools::mixedsort(unique(c(test$i, test$j)))
  test2 <- as.data.frame(test2); names(test2) <- "country"
  test2$hi_m <- test$hi_m[match(test2$country, test$j)]
  test2$hi_x <- test$hi_x[match(test2$country, test$i)]
  test2$V_m <- test$V_m[match(test2$country, test$j)]
  test2$V_x <- test$V_x[match(test2$country, test$i)]
  medM <- medianNA(test2$hi_m)
  medX <- medianNA(test2$hi_x)
  meanM <- meanNA(test2$hi_m)
  meanX <- meanNA(test2$hi_x)
  wmedM <- w.perc(test2$hi_m, test2$V_m, 0.5, na.rm=T)
  wmedX <- w.perc(test2$hi_x, test2$V_x, 0.5, na.rm=T)
  wmeanM <- sumNA(test2$hi_m*test2$V_m)/sumNA(test2$V_m)
  wmeanX <- sumNA(test2$hi_x*test2$V_x)/sumNA(test2$V_x)
  nM <- length(test2$hi_m[is.na(test2$hi_m)==F])
  nX <- length(test2$hi_x[is.na(test2$hi_x)==F])
  wnM <- 1/HHI(test2$V_m)
  wnX <- 1/HHI(test2$V_x)
  
  DF[i,1] <- TovList[i]
  DF[i,2:ncol(DF)] <- cbind(medM,medX,wmedM,wmedX,meanM,meanX,wmeanM,wmeanX,nM,nX,wnM,wnX)
  
  i=i+1
}
##



# 4) Аналитика и регрессии ####
DF$ProdPower <- DF$wmeanM-DF$wmeanX
DF$ProdPower_Med <- DF$wmedM-DF$wmedX


## 4.1) Рис. 2 ####
df_melt <- subset(DF, select=c("TNVED","wmeanM","wmeanX","wnM","wnX"))
df_melt <- reshape2::melt(df_melt)
df_melt$variable <- as.character(df_melt$variable)
# во избежание путаницы меняем обозначения местами (простая 
# концентрация считается напрямую, а двусторонняя - по партнерам):
df_melt$variable[df_melt$variable=="wnM"] <- "wnumX" # это про импорт
df_melt$variable[df_melt$variable=="wnX"] <- "wnumM" # это про экспорт

df_melt$flow <- substr(df_melt$variable,(nchar(df_melt$variable)),nchar(df_melt$variable))
df_melt$var <- gsub("X|M","",df_melt$variable)
df_melt$value[df_melt$var=="wnum"] <- 1/df_melt$value[df_melt$var=="wnum"]

flow1 <- "односторонние потоки"
flow2 <- "двусторонние потоки"
df_cast <- reshape2::dcast(df_melt, TNVED+var~flow, sum, value.var="value")
df_cast$var[df_cast$var=="wmean"] <- flow2
df_cast$var[df_cast$var=="wnum"] <- flow1
df_cast$var <- factor(df_cast$var, c(flow1, flow2))

m2 <- median(df_cast$M[df_cast$var==flow2])
m1 <- median(df_cast$M[df_cast$var==flow1])
x2 <- median(df_cast$X[df_cast$var==flow2])
x1 <- median(df_cast$X[df_cast$var==flow1])

g1 <- ggplot(df_cast) + geom_density(aes(X, col=var, linewidth=var)) + # weight=V
  scale_color_manual(values=c("grey50","black")) + scale_linewidth_manual(values=c(1,0.75)) + 
  xlab("Индекс Херфиндаля-Хиршмена") + ylab("Плотность распределения") + theme_GOST_base() + 
  theme(legend.position="none") + 
  geom_vline(xintercept=x1, col="grey50", linetype="longdash") + 
  geom_vline(xintercept=x2, col="black", linetype="longdash") + 
  ggtitle("Концентрация покупателей") + scale_x_continuous(limits=c(0,1), breaks=seq(0,1,0.2))
g2 <- ggplot(df_cast) + geom_density(aes(M, col=var, linewidth=var)) + # weight=V
  scale_color_manual(values=c("grey50","black")) + scale_linewidth_manual(values=c(1,0.75)) + 
  xlab("Индекс Херфиндаля-Хиршмена") + ylab("Плотность распределения") + theme_GOST_base() + 
  theme(legend.position="bottom", legend.title=element_blank(), legend.text=element_text(size=13)) + 
  geom_vline(xintercept=m1, col="grey50", linetype="longdash") + 
  geom_vline(xintercept=m2, col="black", linetype="longdash") + 
  ggtitle("Концентрация поставщиков") + scale_x_continuous(limits=c(0,1), breaks=seq(0,1,0.2))
png('./pic2.png', 
    width = 6, height = 8, units = 'in', res = 500)
gridExtra::grid.arrange(g1, g2, heights=c(3.7,4.3))
dev.off()
CairoPDF("./pic2.pdf", width = 6.5, height = 8.5)
gridExtra::grid.arrange(g1, g2, heights=c(3.7,4.3))
dev.off()


## 4.2) Табл. 1 ####
tb1 <- as.data.frame(cor(DF[,c("wmedM","wmedX","wmeanM","wmeanX","wnX","wnM","nX","nM")]))
openxlsx::write.xlsx(tb1, "./tab1.xlsx", rowNames=T)


## 4.3) Табл. 2-3 (две регрессии, в 8 вариациях каждая) ####
options(scipen=15)
r=1
while(r<3){
  f1 <- as.formula('scale(wmedM)~scale(wnX)+scale(wnM)+scale(wmedX)')
  f1_ <- as.formula('scale(wmedM)~scale(wnX)+scale(wmedX)')
  f2 <- as.formula('scale(wmeanM)~scale(wnX)+scale(wnM)+scale(wmeanX)')
  f2_ <- as.formula('scale(wmeanM)~scale(wnX)+scale(wmeanX)')
  rf1 <- as.formula('scale(wmedX)~scale(wnM)+scale(wnX)+scale(wmedM)')
  rf1_ <- as.formula('scale(wmedX)~scale(wnM)+scale(wmedM)')
  rf2 <- as.formula('scale(wmeanX)~scale(wnM)+scale(wnX)+scale(wmeanM)')
  rf2_ <- as.formula('scale(wmeanX)~scale(wnM)+scale(wmeanM)')
  
  if(r==1){
    list_f <- list(f1,f1_,f2,f2_)
  } else{
    list_f <- list(rf1,rf1_,rf2,rf2_)
  }
  mods <- list()
  
  j=1
  while(j<length(list_f)+1){
    mods[[j*2-1]] <- lm(formula=list_f[[j]], DF)
    modRLMM1 <- MASS::rlm(list_f[[j]], DF, method="MM") # робастная регрессия
    mods[[j*2]] <- lm(formula=list_f[[j]], DF, weights=modRLMM1$w) # OLS с весами из робастной
    j=j+1
  }
  
  # меняем порядок, чтобы сначала шли 2 простых OLS, а потом 2 робастных
  names(mods) <- glue::glue('col{c(1,3,2,4,5,7,6,8)}')
  mods <- mods[order(names(mods))]
  
  m1 <- "OLS"
  m2 <- "OLS-MM"
  y1 <- ifelse(r==1, "wmedM", "wmedX")
  y2 <- ifelse(r==1, "wmeanM", "wmeanX")
  
  list_coef <- list()
  list_r <- list()
  list_hs <- list()
  list_vif <- list()
  list_p <- list()
  
  j=1
  while(j<length(mods)+1){
    list_coef[[j]] <- round(mods[[j]]$coefficients,3)
    list_r[[j]] <- summary(mods[[j]])$adj.r.squared
    list_hs[[j]] <- as.numeric(lmtest::bptest(mods[[j]])$p.value)
    list_vif[[j]] <- max(car::vif(mods[[j]]))
    list_p[[j]] <- zv(summary(mods[[j]])$coefficients[,4])
    j=j+1
  }
  
  rwNms <- c("Model", unique(c(y1, y2)), unique(Reduce(c, lapply(list_coef, FUN=names))),
             "Adj. R2", "Heterosk. (p-val)", "Multicol. (max VIF)")
  rwNms <- as.data.frame(rwNms)
  j=1
  while(j<length(list_coef)+1){
    rwNms[,j+1] <- list_coef[[j]][match(rwNms$rwNms, names(list_coef[[j]]))]
    j=j+1
  }
  names(rwNms)[2:ncol(rwNms)] <- glue::glue('({1:8})')
  rwNms[(nrow(rwNms)-2),2:ncol(rwNms)] <- round(do.call(c,list_r),3)
  rwNms[(nrow(rwNms)-1),2:ncol(rwNms)] <- round(do.call(c,list_hs),3)
  rwNms[nrow(rwNms),2:ncol(rwNms)] <- round(do.call(c,list_vif),3)
  
  pVals <- rwNms
  j=1
  while(j<length(list_coef)+1){
    pVals[,j+1] <- list_p[[j]][match(pVals$rwNms, names(list_p[[j]]))]
    j=j+1
  }
  
  n <- 2
  while(n < ncol(rwNms)+1){
    rwNms[,n] <- glue::glue('{NAP(rwNms[,{n}])} {NAP(pVals[,{n}])}')
    n=n+1
  }
  
  rwNms[1,c(2,3,6,7)] <- m1
  rwNms[1,c(4,5,8,9)] <- m2
  rwNms[2,2:5] <- "+"
  rwNms[3,6:9] <- "+"
  openxlsx::write.xlsx(rwNms, glue::glue("./tab{r+1}.xlsx"))
  r=r+1
}



# 5) Кластеризация ####
DF_sub <- subset(DF, select=c("TNVED","wmedM","wmedX","wnM","wnX","ProdPower_Med"))
names(DF_sub) <- c("Product","HI2","HJ2","HJ1","HI1","RP")
# для задачи кластеризации показатели эффективного числа 
# экспортеров и импортеров переводятся в индексы концентрации
DF_sub$HJ1 <- 1/DF_sub$HJ1
DF_sub$HI1 <- 1/DF_sub$HI1
DF_sub$SI <- DF_sub$HI2-DF_sub$HI1
DF_sub$SJ <- DF_sub$HJ2-DF_sub$HJ1

results4 <- prcomp(DF_sub[, c(2:5)], scale = TRUE )
results5 <- prcomp(DF_sub[, c(2:6)], scale = TRUE )
results7 <- prcomp(DF_sub[, c(2:8)], scale = TRUE )

prc4 <- round((results4$sdev^2/sum(results4$sdev^2)*100)[1:2],0)
prc5 <- round((results5$sdev^2/sum(results5$sdev^2)*100)[1:2],0)
prc7 <- round((results7$sdev^2/sum(results7$sdev^2)*100)[1:2],0)


## 5.1) Рис. 3 ####
pic_clock4 <- factoextra::fviz_pca_var(
  results4, col.var = "black", alpha.var=0.33, col.circle = alpha("grey70",0.5), repel = TRUE, title="4 переменные", 
  xlab=glue::glue("Первая компонента ({prc4[1]}%)"), ylab=glue::glue("Вторая компонента ({prc4[2]}%)"),
  labelsize = 3) + theme_GOST_base()
pic_clock5 <- factoextra::fviz_pca_var(
  results5, col.var = "black", alpha.var=0.33, col.circle = alpha("grey70",0.5), repel = TRUE, title="5 переменных", 
  xlab=glue::glue("Первая компонента ({prc5[1]}%)"), ylab=glue::glue("Вторая компонента ({prc5[2]}%)"),
  labelsize = 3) + theme_GOST_base()
pic_clock7 <- factoextra::fviz_pca_var(
  results7, col.var = "black", alpha.var=0.33, col.circle = alpha("grey70",0.5), repel = TRUE, title="7 переменных", 
  xlab=glue::glue("Первая компонента ({prc7[1]}%)"), ylab=glue::glue("Вторая компонента ({prc7[2]}%)"),
  labelsize = 3) + theme_GOST_base()

png('./pic3.png', 
    width = 9, height = 3, units = 'in', res = 500)
gridExtra::grid.arrange(pic_clock4, pic_clock5, pic_clock7, ncol=3)
dev.off()
CairoPDF("./pic3.pdf", width = 9.5, height = 3.5)
gridExtra::grid.arrange(pic_clock4, pic_clock5, pic_clock7, ncol=3)
dev.off()


## 5.2) Рис. 4 (оптимальное число кластеров) ####
# (расчет этого блока занимает несколько минут)
LMS <- c(1.05,1.3) # для установки лимитов графиков по оси Y
k_g4 <- fviz_nbclust(results4$x, kmeans, method = "gap_stat", linecolor="black") + scale_y_continuous(limits=LMS) + 
  ggtitle("k-средних, 4 переменные") + xlab("Число кластеров") + ylab("Статистика разрыва") +
  theme_GOST_base()
k_g5 <- fviz_nbclust(results5$x, kmeans, method = "gap_stat", linecolor="black") + scale_y_continuous(limits=LMS) + 
  ggtitle("k-средних, 5 переменных") + xlab("Число кластеров") + ylab("Статистика разрыва") +
  theme_GOST_base()
k_g7 <- fviz_nbclust(results7$x, kmeans, method = "gap_stat", linecolor="black") + scale_y_continuous(limits=LMS) + 
  ggtitle("k-средних, 7 переменных") + xlab("Число кластеров") + ylab("Статистика разрыва") +
  theme_GOST_base()
h_g4 <- fviz_nbclust(results4$x, hcut, method = "gap_stat", linecolor="black") + scale_y_continuous(limits=LMS) + 
  ggtitle("иерархия, 4 переменные") + xlab("Число кластеров") + ylab("Статистика разрыва") +
  theme_GOST_base()
h_g5 <- fviz_nbclust(results5$x, hcut, method = "gap_stat", linecolor="black") + scale_y_continuous(limits=LMS) + 
  ggtitle("иерархия, 5 переменных") + xlab("Число кластеров") + ylab("Статистика разрыва") +
  theme_GOST_base()
h_g7 <- fviz_nbclust(results7$x, hcut, method = "gap_stat", linecolor="black") + scale_y_continuous(limits=LMS) + 
  ggtitle("иерархия, 7 переменных") + xlab("Число кластеров") + ylab("Статистика разрыва") +
  theme_GOST_base()
png('./pic4.png', 
    width = 10, height = 6, units = 'in', res = 500)
gridExtra::grid.arrange(k_g4, k_g5, k_g7, h_g4, h_g5, h_g7, ncol=3)
dev.off()
CairoPDF("./pic4.pdf", width = 10.5, height = 6.5)
gridExtra::grid.arrange(k_g4, k_g5, k_g7, h_g4, h_g5, h_g7, ncol=3)
dev.off()


## 5.3) Табл. 4 (собственно кластеризация) ####
set.seed(13)
clust3 <- kmeans(scale(DF_sub[, 2:8]), 3, iter.max = 100, nstart = 25, algorithm="Hartigan-Wong")
clust3$cluster <- as.factor(clust3$cluster)
DF_sub$clust3 <- as.character(clust3$cluster)
DF_sub$PC1 <- as.data.frame(results7$x)$PC1
DF_sub$PC2 <- as.data.frame(results7$x)$PC2
CLR <- DF_sub %>% group_by(clust3) %>% summarise(
  Kol = sum(as.numeric(clust3)/as.numeric(clust3)), HI2 = median(HI2), HJ2 = median(HJ2), HI1 = median(HI1), HJ1 = median(HJ1),
  RP = median(RP), SI = median(SI), SJ = median(SJ)) %>% as.data.frame()
CLR[,3:ncol(CLR)] <- round(CLR[,3:ncol(CLR)],3)
CLR <- as.data.frame(t(CLR))
MEDS <- DF_sub %>% group_by("всего") %>% summarise(
  Kol = sum(as.numeric(clust3)/as.numeric(clust3)), HI2 = median(HI2), HJ2 = median(HJ2), HI1 = median(HI1), HJ1 = median(HJ1),
  RP = median(RP), SI = median(SI), SJ = median(SJ)) %>% as.data.frame()
MEDS[,3:ncol(MEDS)] <- round(MEDS[,3:ncol(MEDS)],3)
MEDS <- as.data.frame(t(MEDS))
CLR$V4 <- MEDS$V1
openxlsx::write.xlsx(CLR,"./tab4.xlsx", rowNames=T, colNames=F)


## 5.4) Рис. 5 ####
b1 <- ggplot(DF_sub) + geom_text(aes(PC1, PC2, label=clust3, col=clust3), size=2) + theme_GOST_base() +
  xlab(glue::glue("Первая компонента ({prc7[1]}%)")) + ylab(glue::glue("Вторая компонента \n({prc7[2]}%)")) +
  theme(legend.position="none") + scale_color_manual(values=c("grey60","grey40","black","red")) +
  ggtitle("(А)") + scale_y_continuous(breaks=seq(-6,6,2))
b2 <- ggplot(DF_sub) + geom_text(aes(HI2, HJ2, label=clust3, col=clust3), size=2) + theme_GOST_base() +
  xlab(glue::glue("Концентрация поставщиков (HI2)")) + ylab(glue::glue("Концентрация покупателей \n(HJ2)")) +
  theme(legend.position="none") + scale_color_manual(values=c("grey60","grey40","black","red")) +
  ggtitle("(Б)") + scale_x_continuous(breaks=seq(0,1,0.2)) + scale_y_continuous(breaks=seq(0,1,0.2))
b3 <- ggplot(DF_sub) + geom_text(aes(SJ, RP, label=clust3, col=clust3), size=2) + theme_GOST_base() +
  xlab(glue::glue("Сегрегация рынков по покупателям (SJ)")) + 
  ylab(glue::glue("Относительная рыночная \nвласть поставщиков (RP)")) +
  theme(legend.position="none") + scale_color_manual(values=c("grey60","grey40","black","red")) +
  ggtitle("(В)") + scale_x_continuous(breaks=seq(0,1,0.2)) + scale_y_continuous(breaks=seq(-1,1,0.2))
b4 <- ggplot(DF_sub) + geom_text(aes(SI, RP, label=clust3, col=clust3), size=2) + theme_GOST_base() +
  xlab(glue::glue("Сегрегация рынков по поставщикам (SI)")) + 
  ylab(glue::glue("Относительная рыночная \nвласть поставщиков (RP)")) +
  theme(legend.position="none") + scale_color_manual(values=c("grey60","grey40","black","red")) +
  ggtitle("(Г)") + scale_x_continuous(breaks=seq(0,1,0.2)) + scale_y_continuous(breaks=seq(-1,1,0.2))

# комбинированный график из 4 панелей с разделительными линиями между панелями
line_h <- ggplot() + geom_hline(yintercept=1, linetype="solid", linewidth=0.75, col=alpha("black",0.5)) + 
  theme_void()
line_v <- ggplot() + geom_vline(xintercept=1, linetype="solid", linewidth=0.75, col=alpha("black",0.5)) + 
  theme_void() + theme(plot.margin = unit(c(0.25, 0, 0.025, 0), "inches"))
line_no <- ggplot() + theme_void()
part1 <- gridExtra::grid.arrange(b1,line_v,b2, ncol=3, widths=c(8,0.1,8))
part2 <- gridExtra::grid.arrange(b3,line_no,b4, ncol=3, widths=c(8,0.1,8))

png('./pic5.png', 
    width = 10, height = 7, units = 'in', res = 500)
gridExtra::grid.arrange(part1,line_h,part2, heights=c(8,0.1,8))
dev.off()
CairoPDF("./pic5.pdf", width = 10, height = 7)
gridExtra::grid.arrange(part1,line_h,part2, heights=c(8,0.1,8))
dev.off()



# Приложение (состав кластеров и их вклад в торговлю) ####
options(scipen=7)
DF_sub$V <- TOT$v[match(DF_sub$Product, TOT$c)]

ann <- as.data.frame(c(1,2,3))
names(ann) <- "cluster"
j=1
while(j < nrow(ann)+1){
  ann$share[j] <- round(sum(DF_sub$V[DF_sub$clust3==j])/sum(DF_sub$V)*100,1)
  ann$prod_no[j] <- length(DF_sub$Product[DF_sub$clust3==j]) 
  ann$products[j] <- paste(DF_sub$Product[DF_sub$clust3==j], collapse=', ') 
  j=j+1
}
ann$share[1] <- 100-sum(ann$share[2:3])
openxlsx::write.xlsx(ann, "./annex.xlsx")
###
### КОНЕЦ КОДА ###