# <encoding: UTF-8> #

# КОД ДЛЯ ВОСПРОИЗВОДСТВА РЕЗУЛЬТАТОВ СТАТЬИ:
# Гнидченко А.А. (2024). Сдвиги в торговых взаимосвязях США и Китая 
# со странами-партнерами: что изменилось за «пятилетку бурных перемен».
# Современная мировая экономика, Том 2, №2, 2024.


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

## Пользовательские функции
`%!in%` = Negate(`%in%`)
ZNA <- function(x) { x[x==0] <- NA; x }
NAZ <- function(x) { x[is.na(x)==T] <- 0; x }
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)))
  }
}

## Загрузка и подготовка пакетов
packs <- c("openxlsx","reshape2","lubridate","ggplot2","stlplus",
           "data.table","Cairo","grid","gridExtra","readr")
newpacks <- packs[!(packs %in% installed.packages()[,"Package"])]
if(length(newpacks)) install.packages(newpacks)
lapply(packs, library, character.only=T)[[1]]


### 1) Динамика торговли товарами США ####
## Загрузка данных, скачанных отсюда [VPN]: www.census.gov/foreign-trade/balance/country.xlsx
censTB <- openxlsx::read.xlsx("./data/us/country.xlsx")
censTB <- reshape2::melt(censTB, id.vars=c("year", "CTY_CODE", "CTYNAME"))
censTB <- subset(censTB, value!=0)
censTB$Flow <- substr(censTB$variable,1,1)
censTB$mnth <- substr(censTB$variable,2,4)
censTB <- subset(censTB, mnth!="YR")
censTB$mnth <- match(censTB$mnth, toupper(month.abb))
censTB$DATE <- as.Date(paste(censTB$year,censTB$mnth,1,sep="-"))
censTB$DATE <- lubridate::ceiling_date(censTB$DATE, "month")-1
censTB <- reshape2::dcast(censTB, DATE+CTY_CODE+CTYNAME~Flow, sum, value.var="value")
censTB$B <- censTB$E-censTB$I

## Проставление группировок
censTB$CTY_GRP <- NA
eu_cntrs_UsGoods <- c(
  "Austria","Belgium","Bulgaria","Hungary","Germany","Greece","Denmark","Ireland",
  "Spain","Italy","Cyprus","Latvia","Lithuania","Luxembourg","Malta","Netherlands",
  "Poland","Portugal","Romania","Slovakia","Slovenia","Finland","France","Croatia",
  "Czech Republic","Sweden","Estonia")
censTB$CTY_GRP[censTB$CTYNAME%in%c(eu_cntrs_UsGoods,"Norway","Switzerland",
                                   "United Kingdom","Iceland")] <- "Европа"
censTB$CTY_GRP[censTB$CTYNAME%in%c("Vietnam","Malaysia","Indonesia","Thailand","Singapore",
                                   "Philippines","Cambodia","Brunei","Laos","Burma")] <- "АСЕАН"
censTB$CTY_GRP[censTB$CTYNAME%in%c("Canada","Mexico")] <- "Северная Америка"
censTB$CTY_GRP[censTB$CTYNAME%in%c("China")] <- "Китай"
censTB$CTY_GRP[censTB$CTYNAME%in%c("India")] <- "Индия"
censTB$CTY_GRP[censTB$CTYNAME%in%c("Taiwan","Korea, South","Japan")] <- "Альянс Chip4"
censTB$CTY_GRP[censTB$CTYNAME%in%c("World, Not Seasonally Adjusted")] <- "ВСЕГО"

## Создание подвыборки данных с 2015 года
censTB_sub <- subset(censTB, DATE>"2014-12-31" & is.na(CTY_GRP)==F)
censTB_sub$CTY_GRP <- factor(censTB_sub$CTY_GRP, c(
  "ВСЕГО","Китай","Северная Америка","Европа","АСЕАН","Альянс Chip4","Индия"))

## Обработка данных по импорту (в т.ч. расчет тренда)
censTB_cast_m <- reshape2::dcast(censTB_sub, DATE~CTY_GRP, sum, value.var="I")
censTB_cast_m$`Прочие` <- censTB_cast_m$ВСЕГО-apply(censTB_cast_m[,3:ncol(censTB_cast_m)],1,sum)
censTB_melt_m <- reshape2::melt(censTB_cast_m, id.vars=c("DATE"))
names(censTB_melt_m) <- c("DATE","CNTR","VAL")
setDT(censTB_melt_m)[, VAL_SEAS:= stlplus::stlplus(VAL, n.p=12, s.window = 7)$data$trend, by=.(CNTR)]

## Обработка данных по экспорту (в т.ч. расчет тренда)
censTB_cast_x <- reshape2::dcast(censTB_sub, DATE~CTY_GRP, sum, value.var="E")
censTB_cast_x$`Прочие` <- censTB_cast_x$ВСЕГО-apply(censTB_cast_x[,3:ncol(censTB_cast_x)],1,sum)
censTB_melt_x <- reshape2::melt(censTB_cast_x, id.vars=c("DATE"))
names(censTB_melt_x) <- c("DATE","CNTR","VAL")
setDT(censTB_melt_x)[, VAL_SEAS:= stlplus::stlplus(VAL, n.p=12, s.window = 7)$data$trend, by=.(CNTR)]

## Оценка торгового баланса
censTB_melt_m$FLOW <- "импорт"
censTB_melt_x$FLOW <- "экспорт"
censTB_melt_b <- censTB_melt_x
censTB_melt_b$VAL <- NA
censTB_melt_b$VAL_SEAS <- censTB_melt_x$VAL_SEAS-censTB_melt_m$VAL_SEAS
censTB_melt_b$FLOW <- "баланс"
censTB_melt_t <- rbind(censTB_melt_m, censTB_melt_x, censTB_melt_b)

## Построение графика
rect <- data.frame(xmin=as.Date('2018-08-31'), xmax=as.Date('2020-01-31'), ymin=-Inf, ymax=Inf)
fig1 <- ggplot(censTB_melt_t) + 
  geom_line(aes(DATE, VAL/1000, group=interaction(CNTR,FLOW), col=FLOW), alpha = 0.25) +
  geom_line(aes(DATE, VAL_SEAS/1000, group=interaction(CNTR,FLOW), col=FLOW, linewidth=FLOW)) +
  scale_colour_manual(values=c("blue", "darkred","darkgreen")) + 
  scale_linewidth_manual(values=c(1, 0.5, 0.5)) + 
  theme_minimal(base_size=15) + ylab("млрд долл.") + xlab(NULL) + 
  scale_x_date(date_labels="%y", breaks="2 years", expand=c(0.07,0)) +
  facet_wrap(~CNTR, ncol=4, scales="free_y") + 
  geom_hline(yintercept=0, col="darkgrey") + 
  geom_vline(xintercept=as.Date("2022-02-28"), col="red", linetype="dashed", linewidth=0.5, alpha=0.5) +
  geom_rect(data = rect, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax), alpha = 0.2) +
  scale_y_continuous(expand=c(0.1,0)) +
  labs(caption=paste0("Примечание: красный пунктир — начало СВО; ",
                      "серая заливка — активная фаза торговой войны США и Китая")) +
  #ggtitle("Динамика торговли США (01/2015 — 05/2024)") + 
  theme(plot.caption = element_text(size=11, face = "italic"), plot.title = element_text(hjust=0.5),
        axis.text.y = element_text(size=11), legend.box.margin=margin(-10,-10,-10,-10),
        legend.position="bottom", legend.title=element_blank())
fig1

png(paste0('./figs/fig1.png'), 
    width = 8.5, height = 6, units = 'in', res = 500)
fig1
dev.off()

CairoPDF("./figs/fig1.pdf", width = 10.2, height = 7.2)
fig1
dev.off()
###


### 2) Динамика торговли товарами Китая ####
## Загрузка данных, скачанных отсюда: stats.customs.gov.cn/indexEn
FLS_LIST_X <- list.files("./data/chn", pattern="exW", full.names=T)
chn_exW <- rbindlist(lapply(
  FLS_LIST_X[], function(k){read.csv(k, as.is=T, sep=",", na.strings = "")}))
chn_exW$Flow <- "Export"
FLS_LIST_M <- list.files("./data/chn", pattern="imW", full.names=T)
chn_imW <- rbindlist(lapply(
  FLS_LIST_M[], function(k){read.csv(k, as.is=T, sep=",", na.strings = "")}))
chn_imW$Flow <- "Import"

## Обработка данных
chn_trW <- rbind(chn_exW, chn_imW)
chn_trW$US.dollar <- as.numeric(gsub(",","",chn_trW$US.dollar))/1000000
chn_trW$Date <- as.Date(paste(substr(
  chn_trW$Date.of.data,1,4), substr(chn_trW$Date.of.data,5,6), "01", sep="-"))
chn_trW$Date <- lubridate::ceiling_date(chn_trW$Date, "month")-1

## Расчет баланса
chn_trW_cast <- reshape2::dcast(
  chn_trW, Date+Trading.partner.code+Trading.partner~Flow, sum, value.var="US.dollar")
chn_trW_cast$Balance <- chn_trW_cast$Export-chn_trW_cast$Import
names(chn_trW_cast) <- c("DATE","CTY_CODE","CTY_NAME","X","I","B")

## Корректировка некоторых названий стран
chn_trW_cast$CTY_NAME[chn_trW_cast$CTY_NAME=="T\xa8\xb9rkiye"] <- "Turkiye"
chn_trW_cast$CTY_NAME <- iconv(chn_trW_cast$CTY_NAME, from="GB2312", to="UTF-8")
chn_trW_cast$CTY_NAME[chn_trW_cast$CTY_NAME=="Cura?ao"] <- "Curacao"
chn_trW_cast$CTY_NAME <- gsub("\\?","A",chn_trW_cast$CTY_NAME)
chn_trW_cast$CTY_NAME <- gsub("[\r\n]","",chn_trW_cast$CTY_NAME)
chn_trW_cast$CTY_CODE <- as.character(chn_trW_cast$CTY_CODE)

## Расчет агрегата торговли со всеми странами
chn_tb_agr <- aggregate(cbind(X,I,B)~DATE, chn_trW_cast, sum)
chn_tb_agr$CTY_CODE <- "000"
chn_tb_agr$CTY_NAME <- "World"
chn_tb_agr <- subset(chn_tb_agr, select=names(chn_trW_cast))
chn_trW_cast <- rbind(chn_trW_cast, chn_tb_agr)

## Проставление группировок
chn_trW_cast$CTY_GRP <- NA
eu_cntrs_CnGoods <- c(
  "Austria","Belgium","Bulgaria","Hungary","Germany","Greece","Denmark","Ireland",
  "Spain","Italy","Cyprus","Latvia","Lithuania","Luxembourg","Malta","Netherlands",
  "Poland","Portugal","Romania","Slovakia","Slovenia","Finland","France","Croatia",
  "Czechia(Czech Republic)","Sweden","Estonia")
chn_trW_cast$CTY_GRP[chn_trW_cast$CTY_NAME%in%c(eu_cntrs_CnGoods,"Norway","Switzerland",
                                   "United Kingdom","Iceland")] <- "Европа"
chn_trW_cast$CTY_GRP[chn_trW_cast$CTY_NAME%in%c("World")] <- "ВСЕГО"
chn_trW_cast$CTY_GRP[chn_trW_cast$CTY_NAME%in%c("United States")] <- "США"
chn_trW_cast$CTY_GRP[chn_trW_cast$CTY_NAME%in%c("Canada","Mexico")] <- "Северная Америка"
chn_trW_cast$CTY_GRP[chn_trW_cast$CTY_NAME%in%c("India")] <- "Индия"
chn_trW_cast$CTY_GRP[chn_trW_cast$CTY_NAME%in%c("Japan","Republic of Korea","Taiwan,China")] <- "Альянс Chip4"
chn_trW_cast$CTY_GRP[chn_trW_cast$CTY_NAME%in%c("Viet Nam","Malaysia","Indonesia","Thailand","Singapore",
                                                "Philippines","Cambodia","Brunei","Lao","Myanmar")] <- "АСЕАН"

##  Корректировака на реимпорт (Гонконг)
reimp_look <- as.data.frame(cbind(
  chn_trW_cast$DATE[chn_trW_cast$CTY_NAME=="World"],
  chn_trW_cast$I[chn_trW_cast$CTY_NAME=="World"],
  chn_trW_cast$X[chn_trW_cast$CTY_NAME=="World"],
  chn_trW_cast$I[chn_trW_cast$CTY_NAME=="China"]))
reimp_look$V1 <- as.Date(reimp_look$V1, origin="1970-01-01")
names(reimp_look) <- c("date","imp","exp","reimp")

chn_trW_cast$REIMP <- chn_trW_cast$I[match(
  interaction(chn_trW_cast$DATE, "China"),
  interaction(chn_trW_cast$DATE, chn_trW_cast$CTY_NAME))]
chn_trW_cast$X[chn_trW_cast$CTY_NAME=="Hong Kong,China"] <- 
  chn_trW_cast$X[chn_trW_cast$CTY_NAME=="Hong Kong,China"]-
  chn_trW_cast$REIMP[chn_trW_cast$CTY_NAME=="Hong Kong,China"]
chn_trW_cast$I[chn_trW_cast$CTY_NAME=="China"] <- 
  chn_trW_cast$I[chn_trW_cast$CTY_NAME=="China"]-
  chn_trW_cast$REIMP[chn_trW_cast$CTY_NAME=="Hong Kong,China"]
chn_trW_cast$X[chn_trW_cast$CTY_NAME=="World"] <- 
  chn_trW_cast$X[chn_trW_cast$CTY_NAME=="World"]-
  chn_trW_cast$REIMP[chn_trW_cast$CTY_NAME=="Hong Kong,China"]
chn_trW_cast$I[chn_trW_cast$CTY_NAME=="World"] <- 
  chn_trW_cast$I[chn_trW_cast$CTY_NAME=="World"]-
  chn_trW_cast$REIMP[chn_trW_cast$CTY_NAME=="Hong Kong,China"]
chn_trW_cast$B <- chn_trW_cast$X-chn_trW_cast$I

## Подготовка к построению графика
chn_trW_sub <- subset(chn_trW_cast, is.na(CTY_GRP)==F)
chn_trW_sub$CTY_GRP <- factor(chn_trW_sub$CTY_GRP, c(
  "ВСЕГО","США","Северная Америка","Европа","АСЕАН","Альянс Chip4","Индия"))

## Обработка данных по импорту (в т.ч. построение тренда)
chn_trW_cast2_m <- reshape2::dcast(chn_trW_sub, DATE~CTY_GRP, sum, value.var="I")
chn_trW_cast2_m$`Прочие` <- chn_trW_cast2_m$ВСЕГО - 
  apply(chn_trW_cast2_m[,3:ncol(chn_trW_cast2_m)],1,sum)
chn_trW_melt_m <- reshape2::melt(chn_trW_cast2_m, id.vars=c("DATE"))
names(chn_trW_melt_m) <- c("DATE","CNTR","VAL")
setDT(chn_trW_melt_m)[, VAL_SEAS:= stlplus::stlplus(
  VAL, n.p=12, s.window = 7)$data$trend, by=.(CNTR)]

## Обработка данных по экспорту (в т.ч. построение тренда)
chn_trW_cast2_x <- reshape2::dcast(chn_trW_sub, DATE~CTY_GRP, sum, value.var="X")
chn_trW_cast2_x$`Прочие` <- chn_trW_cast2_x$ВСЕГО-apply(chn_trW_cast2_x[,3:ncol(chn_trW_cast2_x)],1,sum)
chn_trW_melt_x <- reshape2::melt(chn_trW_cast2_x, id.vars=c("DATE"))
names(chn_trW_melt_x) <- c("DATE","CNTR","VAL")
setDT(chn_trW_melt_x)[, VAL_SEAS:= stlplus::stlplus(VAL, n.p=12, s.window = 7)$data$trend, by=.(CNTR)]

## Оценка торгового баланса
chn_trW_melt_m$FLOW <- "импорт"
chn_trW_melt_x$FLOW <- "экспорт"
chn_trW_melt_b <- chn_trW_melt_x
chn_trW_melt_b$VAL <- NA
chn_trW_melt_b$VAL_SEAS <- chn_trW_melt_x$VAL_SEAS-chn_trW_melt_m$VAL_SEAS
chn_trW_melt_b$FLOW <- "баланс"
chn_trW_melt_t <- rbind(chn_trW_melt_m, chn_trW_melt_x, chn_trW_melt_b)

## Построение графика
rect <- data.frame(xmin=as.Date('2018-08-31'), xmax=as.Date('2020-01-31'), ymin=-Inf, ymax=Inf)
fig2 <- ggplot(chn_trW_melt_t) + 
  geom_line(aes(DATE, VAL/1000, group=interaction(CNTR,FLOW), col=FLOW), alpha = 0.25) +
  geom_line(aes(DATE, VAL_SEAS/1000, group=interaction(CNTR,FLOW), col=FLOW, linewidth=FLOW)) +
  scale_colour_manual(values=c("blue", "darkred","darkgreen")) + 
  scale_linewidth_manual(values=c(1, 0.5, 0.5)) + 
  theme_minimal(base_size=15) + 
  scale_x_date(date_labels="%y", breaks="2 years", expand=c(0.07,0)) +
  ylab("млрд долл.") + xlab(NULL) + facet_wrap(~CNTR, ncol=4, scales="free_y") + 
  geom_hline(yintercept=0, col="darkgrey") + 
  geom_vline(xintercept=as.Date("2022-02-28"), col="red", linetype="dashed", linewidth=0.5, alpha=0.5) +
  geom_rect(data = rect, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax), alpha = 0.2) +
  scale_y_continuous(expand=c(0.1,0)) +
  labs(caption=paste0("Примечание: красный пунктир — начало СВО; ",
                      "серая заливка — активная фаза торговой войны США и Китая")) +
  #ggtitle("Динамика торговли Китая (01/2015 — 06/2024)") + 
  theme(plot.caption = element_text(size=11, face = "italic"), plot.title = element_text(hjust=0.5),
        axis.text.y = element_text(size=11), legend.box.margin=margin(-10,-10,-10,-10),
        legend.position="bottom", legend.title=element_blank())
fig2

png(paste0('./figs/fig2.png'), 
    width = 8.5, height = 6, units = 'in', res = 500)
fig2
dev.off()

CairoPDF("./figs/fig2.pdf", width = 10.2, height = 7.2)
fig2
dev.off()
###


### 3) Доля Китая и США в мировой торговле ####
## Загрузка данных, скачанных отсюда [VPN]:
# https://www.cpb.nl/sites/default/files/omnidownload/CPB-World-Trade-Monitor-June-2024.xlsx
cpb <- openxlsx::read.xlsx("./data/cpb/CPB-World-Trade-Monitor-June-2024.xlsx", sheet="trade_out")

## Первичная обработка данных
per_nms <- as.character(cpb[2,4:ncol(cpb)])
cpb[,4:ncol(cpb)] <- apply(cpb[,4:ncol(cpb)],2,as.numeric)
cpb[4:nrow(cpb),4:ncol(cpb)] <- cpb[4:nrow(cpb),4:ncol(cpb)]/100

v10 <- cpb[row.names(cpb[cpb[,2]%like%"qnmi_sn",]),c(1,2,3)]
v10[,3] <- as.numeric(v10[,3]) # объем торговли в базовом 2010 году
v_index <- cpb[row.names(cpb[cpb[,2]%like%"qnmi_sn",]),c(1,2,4:ncol(cpb))]
p_index <- cpb[row.names(cpb[cpb[,2]%like%"pdmi_sn",]),c(1,2,4:ncol(cpb))]

## Расчет стоимостных объемов торговли в постоянных и текущих ценах
vols <- v_index # постоянные цены
vols[,3:ncol(vols)] <- v_index[,3:ncol(v_index)]*v10[,3]
names(vols) <- c("indicator","code",per_nms)
vals <- v_index # текущие цены
vals[,3:ncol(vals)] <- v_index[,3:ncol(v_index)]*p_index[,3:ncol(p_index)]*v10[,3]
names(vals) <- c("indicator","code",per_nms)
x <- vals[row.names(vals[vals[,2]%like%"xgz",]),]
m <- vals[row.names(vals[vals[,2]%like%"mgz",]),]

## Расчет доли США и Китая в мировом экспорте и импорте
shrs <- rbind(x[x[,2]%like%"us",3:ncol(x)]/x[1,3:ncol(x)],
              x[x[,2]%like%"cn",3:ncol(x)]/x[1,3:ncol(x)],
              m[m[,2]%like%"us",3:ncol(m)]/m[1,3:ncol(m)],
              m[m[,2]%like%"cn",3:ncol(m)]/m[1,3:ncol(m)])
shrs$cntr <- c("США","Китай","США","Китай")
shrs$flow <- c("доля в экспорте","доля в экспорте","доля в импорте","доля в импорте")

## Обработка таблицы перед построением графика
shrs <- reshape2::melt(shrs, id.vars=c("cntr","flow"))
names(shrs)[3] <- "date"
shrs$date <- as.Date(glue::glue("{gsub('m','-',shrs$date)}-01"))
shrs <- subset(shrs, shrs$date>"2016-12-01")

## Построение графика
fig3 <- ggplot(shrs) + geom_line(aes(date, value*100, col=flow)) + facet_wrap(~cntr) +
  theme_minimal(base_size=16) + xlab(NULL) + ylab(NULL) + 
  theme(legend.position="bottom", legend.title=element_blank(),
        legend.box.margin=margin(-10,-10,-10,-10)) +
  scale_color_manual(values=c("firebrick","darkgreen")) +
  scale_y_continuous(expand=c(0,2), breaks=seq(0,20,2)) + 
  scale_x_date(date_labels="%y", breaks="2 years", expand=c(0.07,0))
fig3

png(paste0('./figs/fig3.png'), 
    width = 8, height = 3, units = 'in', res = 500)
fig3
dev.off()

CairoPDF("./figs/fig3.pdf", width = 9.6, height = 3.6)
fig3
dev.off()
#


### 4) Роль России в экспорте Китая в прочие страны ####
## Создание выборок данных по экспорту Китая в прочие страны и в Россию
dRest <- chn_trW_melt_x[chn_trW_melt_x$CNTR=="Прочие",]
dRuss <- chn_trW_cast[chn_trW_cast$CTY_NAME=="Russia",1:4]
dRuss$VAL_SEAS <- stlplus::stlplus(dRuss$X, n.p=12, s.window = 7)$data$trend

## Склейка двух подвыборок в одну таблицу
dRest$VAL_RU <- dRuss$X # Россия
dRest$VAL_SEAS_RU <- dRuss$VAL_SEAS
dRest$VAL_krRU <- dRest$VAL-dRest$VAL_RU # кроме России
dRest$VAL_SEAS_krRU <- dRest$VAL_SEAS-dRest$VAL_SEAS_RU
dRest <- subset(dRest, DATE>"2020-11-30") # выборка данных с декабря 2020 г.

## Обработка данных для построения первой части рисунка
d_f4a <- subset(dRest, select=c("DATE","VAL","VAL_SEAS","VAL_krRU","VAL_SEAS_krRU"))
d_f4a <- reshape2::melt(d_f4a, id.vars="DATE")
d_f4a$variable <- as.character(d_f4a$variable)
d_f4a$variable[d_f4a$variable=="VAL_SEAS"] <- "ВСЕГО (тренд)"
d_f4a$variable[d_f4a$variable=="VAL_SEAS_krRU"] <- "кроме России (тренд)"
d_f4a$variable[d_f4a$variable=="VAL"] <- "ВСЕГО (исходные данные)"
d_f4a$variable[d_f4a$variable=="VAL_krRU"] <- "кроме России (исходные данные)"

## Построение первой части рисунка
fig4_1 <- ggplot(d_f4a[d_f4a$DATE>"2020-12-31",]) + 
  geom_line(aes(DATE, value/1000, group=variable, col=variable, linewidth=variable)) + 
  theme_minimal(base_size=20) + xlab(NULL) + ylab(NULL) + 
  scale_color_manual(values=c(alpha("red",0.33),"red",alpha("orange",0.33),"orange")) +
  scale_linewidth_manual(values=c(0.5, 1, 0.5, 1)) + guides(col=guide_legend(nrow=2, byrow=T)) +
  theme(legend.position="bottom", legend.title=element_blank(), legend.box.margin=margin(-10,-10,-10,-10),
        plot.title = element_text(hjust=0.5)) +
  ggtitle("Объем экспорта Китая \nв прочие страны (по тренду)")
fig4_1

## Обработка данных для построения второй части рисунка
d_f4b <- reshape2::dcast(d_f4a[d_f4a$variable%like%"тренд",], DATE~variable, sum, value.var="value")
names(d_f4b)[2:3] <- gsub(" [(]тренд[)]","",names(d_f4b)[2:3])
d_f4b$`ВСЕГО` <- c(NA, d_f4b$`ВСЕГО`[2:nrow(d_f4b)]-d_f4b$`ВСЕГО`[1:(nrow(d_f4b)-1)])
d_f4b$`кроме России` <- c(NA, d_f4b$`кроме России`[2:nrow(d_f4b)]-d_f4b$`кроме России`[1:(nrow(d_f4b)-1)])
d_f4b$`Россия` <- d_f4b$`ВСЕГО`-d_f4b$`кроме России`
d_f4b <- d_f4b[-1,]
d_f4b <- reshape2::melt(d_f4b, id.vars=c("DATE"))
d_f4b$variable <- factor(d_f4b$variable, c("Россия","кроме России","ВСЕГО"))

## Построение второй части рисунка
fig4_2 <- ggplot() + 
  geom_col(data=d_f4b[d_f4b$variable%in%c("Россия","кроме России"),],aes(
    DATE, value/1000, group=variable, fill=variable), position="stack") +
  scale_fill_manual(values=c(alpha("blue",0.2),alpha("orange",0.33))) + 
  geom_point(data=d_f4b[d_f4b$variable%in%c("ВСЕГО"),], aes(DATE, value/1000, col=variable), size=1.5) +
  geom_col(data=d_f4b[d_f4b$variable%in%c("ВСЕГО"),], aes(DATE, value/1000), fill="firebrick", width=0.75) +
  scale_color_manual(values=c("firebrick")) + theme_minimal(base_size=20) + xlab(NULL) + ylab(NULL) + 
  theme(legend.box.margin=margin(-10,-10,-10,-10), legend.position="bottom", legend.title=element_blank(),
        plot.title = element_text(hjust=0.5)) +
  ggtitle("Прирост экспорта Китая \nв прочие страны (по тренду)")
fig4_2

## Построение рисунка в целом
library("gridExtra")
png(paste0('./figs/fig4.png'), 
    width = 14, height = 5, units = 'in', res = 500)
grid.arrange(fig4_1,fig4_2, ncol=2)
dev.off()

CairoPDF("./figs/fig4.pdf", width = 16.8, height = 6)
grid.arrange(fig4_1,fig4_2, ncol=2)
dev.off()
####


## 5) Торговля услугами США ####
## Загрузка данных по мировой торговле услугами с сайта ВТО
if(file.exists("./data/serv/services_annual_dataset.zip")==F){
  url <- "https://stats.wto.org/assets/UserGuide/services_annual_dataset.zip"
  download.file(url, "./data/serv/services_annual_dataset.zip", method="curl", mode="wb")
}

## Получение таблицы из загруженного zip-архива и подвыборка данных с 2015 г.
serv <- read_csv("./data/serv/services_annual_dataset.zip")
serv <- subset(serv, serv$Year>2014)

## Сохранение данных по кодам услуг, не участвующим в основных расчетах
serv_agr <- subset(serv, serv$ProductCode%in%c("S","SA","SJ","SN"))
serv_agr <- subset(serv_agr, serv_agr$ReporterISO3A=="USA")
serv_agr <- subset(serv_agr, serv_agr$Partner!="European Union")
serv_agr$Partner[serv_agr$Partner=="T\xfcrkiye"] <- "Turkiye"

## Подвыборка данных по кодам услуг, участвующим в основных расчетах
serv <- subset(serv, serv$ProductCode%in%c(
  "SA","SB","SC","SD","SE","SF","SG","SH","SI","SJ1","SJ2","SJ3","SK","SL"))
serv <- subset(serv, serv$ReporterISO3A=="USA")
serv <- subset(serv, serv$Partner!="European Union")
serv$Partner[serv$Partner=="T\xfcrkiye"] <- "Turkiye"
gc()

## Установка страновых групп
serv$CTY_GRP <- NA
serv$CTY_GRP[serv$Partner%in%c("Viet Nam","Malaysia","Indonesia","Thailand","Singapore",
                                           "Philippines","Cambodia","Brunei Darussalam",
                                           "Lao People's Democratic Republic","Myanmar")] <- "АСЕАН"
serv$CTY_GRP[serv$Partner%in%c("Canada","Mexico")] <- "Северная Америка"
serv$CTY_GRP[serv$Partner%in%c("China")] <- "Китай"
eu_cntrs <- c("Austria","Belgium","Bulgaria","Hungary","Germany","Greece","Denmark","Ireland",
              "Spain","Italy","KIPR","Latvia","Lithuania","Luxembourg","Malta","Netherlands",
              "Poland","Portugal","Romania","Slovakia","Slovenia","Finland","France","Croatia",
              "Czechia(Czech Republic)","Sweden","Estonia")
serv$CTY_GRP[serv$Partner%in%c(eu_cntrs,"Norway","Switzerland",
                                           "United Kingdom","Iceland")] <- "Европа"
serv$CTY_GRP[serv$Partner%in%c("India")] <- "Индия"
serv$CTY_GRP[serv$Partner%in%c("Chinese Taipei","Korea, Republic of","Japan")] <- "Альянс Chip4"
serv$CTY_GRP[serv$Partner%in%c("World")] <- "ВСЕГО"

## Подвыборка только по странам заданных групп
serv_sub <- subset(serv, is.na(CTY_GRP)==F)
serv_sub$CTY_GRP <- factor(serv_sub$CTY_GRP, c(
  "ВСЕГО","Китай","Северная Америка","Европа","АСЕАН","Альянс Chip4","Индия"))
serv_sub$ProductCode <- factor(serv_sub$ProductCode, sort(unique(serv_sub$ProductCode)))
serv_sub <- as.data.frame(serv_sub)

## Подготовка переменных для цикла
subs <- list(as.character(unique(serv_sub$ProductCode)),
             c("SB","SE","SF","SG","SH","SI","SJ","SK","SL"),
             "SC","SD",c("SF","SG"),"SI","SH","SJ1")
lbls <- list("Услуги, всего","Услуги, кроме транспортных и туристических",
             "Транспортные услуги","Туристические услуги","Страховые и финансовые услуги",
             "Телеком, компьютерные и информационные услуги",
             "Плата за пользование интеллектуальной собственностью",
             "Услуги исследований и разработок")
fig_names <- c("5","6",glue::glue("A{1:6}"))

## Вспомогательный этап: закрытие дыр в исходных данных для США
# (по каждой стране-партнеру раскидываем дыры структурой прошлого 
# года по видам услуг, последовательно двигаясь от 2016 г. к 2023 г.)
lst_f <- list()
flws <- unique(serv_sub$IndicatorCode)
pars <- sort(unique(serv_sub$PartnerCode))
# [начало вложенного цикла]
f <- 1 # торговый поток: импорт/экспорт (trade flow)
while(f<length(flws)+1){
  tst_lst <- list()
  i=1
  while(i<length(pars)+1){
    # выбор страны-партнера, для которой затем производится расчет
    par <- pars[i]
    tst <- subset(serv_sub[,1:(ncol(serv_sub)-1)], PartnerCode==par & IndicatorCode==flws[f])
    tst_agr <- subset(serv_agr, PartnerCode==par & IndicatorCode==flws[f])
    tst <- rbind(tst_agr, tst)
    
    # закрытие дыр на уровне крупных видов услуг
    tst_0 <- tst[tst$ProductCode%!in%c("SJ1","SJ2","SJ3"),]
    tst_0 <- reshape2::dcast(tst_0, PartnerCode+PartnerISO3A+ProductCode~Year, sum, value.var="Value")
    tst_0[,4:ncol(tst_0)] <- ZNA(tst_0[,4:ncol(tst_0)])
    ttl <- tst_0[1,4:ncol(tst_0)]
    dfr <- ttl-apply(tst_0[2:nrow(tst_0),4:ncol(tst_0)],2,sumNA)
    years_full <- names(dfr)
    tst_1 <- tst_0
    n=2 # старт мини-цикла с 2016 г.
    while(n<length(years_full)+1){
      y <- years_full[n]
      y_prv <- as.character(as.numeric(y)-1)
      if(dfr[n]/ttl[n]*100 > 0.1){
        c1 <- tst_0[2:nrow(tst_0),y]
        raspr <- tst_0[1,y]-sumNA(c1)
        flags <- 1-NAZ(c1/c1)
        c2 <- NAZ((tst_0[2:nrow(tst_0),y_prv]*flags)/
                    sumNA((tst_0[2:nrow(tst_0),y_prv]*flags)))
        c3 <- ZNA(NAZ(c1) + round(raspr*c2,0))
        tst_1[2:nrow(tst_0),y] <- c3
      }
      n=n+1
    }
    
    # закрытие дыр внутри группы SJ
    tst_sj <- tst[substr(tst$ProductCode,1,2)=="SJ",]
    tst_sj <- reshape2::dcast(tst_sj, PartnerCode+PartnerISO3A+ProductCode~Year, sum, value.var="Value")
    tst_sj[tst_sj$ProductCode=="SJ",] <- tst_1[tst_1$ProductCode=="SJ",]
    tst_sj[,4:ncol(tst_sj)] <- ZNA(tst_sj[,4:ncol(tst_sj)])
    ttl <- tst_sj[1,4:ncol(tst_sj)]
    dfr <- ttl-apply(tst_sj[2:nrow(tst_sj),4:ncol(tst_sj)],2,sumNA)
    
    tst_sj_1 <- tst_sj
    n=2 # старт мини-цикла с 2016 г.
    while(n<length(years_full)+1){
      y <- years_full[n]
      y_prv <- as.character(as.numeric(y)-1)
      if((dfr[n]/ttl[n]*100 > 0.1) & (is.na(ttl[n])==F)){
        c1 <- tst_sj_1[2:nrow(tst_sj_1),y]
        raspr <- tst_sj_1[1,y]-sumNA(c1)
        flags <- 1-NAZ(c1/c1)
        c2 <- NAZ((tst_sj_1[2:nrow(tst_sj_1),y_prv]*flags)/
                    sumNA((tst_sj_1[2:nrow(tst_sj_1),y_prv]*flags)))
        c3 <- ZNA(NAZ(c1) + round(raspr*c2,0))
        tst_sj_1[2:nrow(tst_sj_1),y] <- c3
      }
      n=n+1
    }
    
    # сборка исправленной таблицы по конкретной стране-партнеру
    # (только в части кодов услуг, участвующих в основных расчетах)
    tst_agr <- tst_1[tst_1$ProductCode=="S",]
    tst_1 <- tst_1[tst_1$ProductCode%!in%c("S","SJ"),]
    tst_1 <- rbind(tst_1, tst_sj_1[2:nrow(tst_sj_1),])
    tst_1 <- tst_1[order(tst_1$ProductCode),]
    tst_melt <- reshape2::melt(tst_1, id.vars=names(tst_1)[1:3])
    names(tst_melt)[4:5] <- c("Year","Value")
    
    # накопление исправленных таблиц в листе (для конкретного торгового потока)
    tst_lst[[i]] <- tst_melt
    i=i+1
  }
  
  # построение исправленной таблицы для всего торгового потока
  tst_df <- do.call(rbind, tst_lst)
  tst_df$IndicatorCode <- flws[f]
  
  # накопление исправленных таблиц в листе (для всех торговых потоков)
  lst_f[[f]] <- tst_df
  f=f+1
}
# [конец вложенного цикла]

## Построение итоговой исправленной таблицы
tst_df <- do.call(rbind, lst_f)

## Обработка итоговой таблицы перед построением графиков
tst_df$Value <- NAZ(tst_df$Value)
tst_df$CTY_GRP <- serv_sub$CTY_GRP[match(tst_df$PartnerCode, serv_sub$PartnerCode)]

## Построение графиков (цикл)
i=1
while(i<length(lbls)+1){
  serv_x <- reshape2::dcast(
    tst_df[tst_df$ProductCode%in%subs[[i]] & tst_df$IndicatorCode=="ITS_CS_AX6",],
    Year~CTY_GRP, sum, value.var="Value")
  serv_x$`Прочие` <- serv_x$ВСЕГО-apply(serv_x[,3:ncol(serv_x)],1,sum)
  serv_x_melt <- reshape2::melt(serv_x, id.vars=c("Year"))
  names(serv_x_melt) <- c("DATE","CNTR","VAL")
  
  serv_m <- reshape2::dcast(
    tst_df[tst_df$ProductCode%in%subs[[i]] & tst_df$IndicatorCode=="ITS_CS_AM6",],
    Year~CTY_GRP, sum, value.var="Value")
  serv_m$`Прочие` <- serv_m$ВСЕГО-apply(serv_m[,3:ncol(serv_m)],1,sum)
  serv_m_melt <- reshape2::melt(serv_m, id.vars=c("Year"))
  names(serv_m_melt) <- c("DATE","CNTR","VAL")
  
  serv_t <- serv_x_melt
  names(serv_t)[3] <- "экспорт"
  serv_t$`импорт` <- serv_m_melt$VAL
  serv_t$`баланс` <- serv_t$`экспорт`-serv_t$`импорт`
  
  serv_t_melt <- reshape2::melt(serv_t, id.vars=c("DATE","CNTR"))
  names(serv_t_melt)[3:4] <- c("FLOW","VAL")
  serv_t_melt$FLOW <- factor(serv_t_melt$FLOW, c("баланс","импорт","экспорт"))
  
  pic_US_serv <- ggplot(serv_t_melt) + 
    theme_minimal(base_size=15) + geom_hline(yintercept=0, col="darkgrey") +  
    geom_line(aes(as.numeric(substr(DATE,3,4)), VAL/1000, group=interaction(CNTR,FLOW), col=FLOW, linewidth=FLOW)) +
    scale_colour_manual(values=c("blue", "darkred","darkgreen")) + 
    scale_linewidth_manual(values=c(1, 0.5, 0.5)) + 
    ylab("млрд долл.") + xlab(NULL) + 
    scale_x_continuous(expand=c(0.07,0), breaks=seq(15,23,2)) +
    scale_y_continuous(expand=c(0.1,0)) + 
    theme(legend.position="bottom", legend.title=element_blank()) + 
    facet_wrap(~CNTR, ncol=4, scales="free_y") + #ggtitle(lbls[[i]]) + 
    theme(plot.caption = element_text(size=11, face = "italic"), 
          plot.title = element_text(hjust=0.5),
          axis.text.y = element_text(size=11), 
          legend.box.margin=margin(-10,-10,-10,-10))
  pic_US_serv
  
  png(glue::glue('./figs/fig{fig_names[[i]]}.png'), 
      width = 8.5, height = 6, units = 'in', res = 500)
  print(pic_US_serv)
  dev.off()
  
  CairoPDF(glue::glue('./figs/fig{fig_names[[i]]}.pdf'),
           width = 10.2, height = 7.2)
  print(pic_US_serv)
  dev.off()
  
  i=i+1
}
### конец кода