analysis/01_baselines.R

In okrebs/DeGVC: Decoupling Global Value Chains companion package

# cf_baselines.R ###############################################################
#                                                                              #
# This script calculates baseline scenarios:                                   #
# 1. without global value chains, i.e. if only intermediate trade is prevented #
#    through infinite intermediate good trade barriers                         #
# 2. with increased intermediate trade barriers                                #
# 3. with overall increased trade barriers                                     #
#                                                                              #
################################################################################

# Libraries ====================================================================

library(tidyverse)
library(haven)
library(DeGVC)
library(here)

# Set up =======================================================================

# Data =========================================================================

# Part of DeGVC pacakge
data("trade_elasticities")
data("intersec_elas")

# Needs to be created by the use, file size too large for github
if(!file.exists(here("data", "wiot.rda"))) {
  stop("First run the script(s) in the folder data-raw!")
} else {
  data(wiot)
}

# Derive simulation inputs -----------------------------------------------------

wiot <- mutate(wiot, simulation_data = map(iot, convert_io))

# Define shocks ================================================================

locations <- wiot$simulation_data[[1]]$location_id # always the same

J <- length(locations)
S <- length(wiot$simulation_data[[1]]$sector_id) # always the same

CHN_id <- match("CHN", locations)
US_id  <- match("USA", locations)
EU_id <- match(c("AUT", "BEL", "BGR", "CYP", "CZE", "DEU", "DNK", "ESP", "EST",
                 "FIN", "FRA", "GBR", "GRC", "HRV", "HUN", "IRL", "ITA", "LTU",
                 "LUX", "LVA", "MLT", "NLD", "POL", "PRT", "ROU", "SVK", "SVN",
                 "SWE"), locations)

shock <- list()

# no change --------------------------------------------------------------------


shock$unchanged <- list(
  T_hat = matrix(1, nrow = J, ncol = S),
  tau_hat = rep(1, J * J * S * (S + 1)),
  delta_hat = matrix(1, nrow = J, ncol = S),
  varkappa_hat = rep(1, J)
)

# gvc --------------------------------------------------------------------------

# all shocks are initially defined as 1% shocks with the magnitude changed later
tau_hat <- matrix(1.01, nrow = J, ncol = J)
diag(tau_hat) <- 1

shock$gvc <- shock$unchanged
shock$gvc$tau_hat <- c(rep(tau_hat, S * S),
                       rep(1, J * J * S)) # final goods barriers unchanged

# final ------------------------------------------------------------------------

shock$final <- shock$unchanged
shock$final$tau_hat <- c(rep(1, J * J * S * S), # intermediates unchanged
                         rep(tau_hat, S))

# total  -----------------------------------------------------------------------

shock$total <- shock$unchanged
shock$total$tau_hat <- rep(tau_hat, S * (S + 1))

# gvc except EU ----------------------------------------------------------------

tau_hat <- matrix(1.01, nrow = J, ncol = J)
diag(tau_hat) <- 1
tau_hat[EU_id, EU_id] <- 1 # not between EU members

shock$gvc_except_eu <- shock$unchanged
shock$gvc_except_eu$tau_hat <- c(rep(as.vector(tau_hat), S * S),
                                 rep(1, J * J * S)) # final unchanged

# final except EU --------------------------------------------------------------

shock$final_except_eu <- shock$unchanged
shock$final_except_eu$tau_hat <- c(rep(1, J * J * S * S), # interm. unchanged
                                   rep(tau_hat, S))

# total except EU --------------------------------------------------------------

shock$total_except_eu <- shock$unchanged
shock$total_except_eu$tau_hat <- rep(tau_hat, S * (S + 1))

# bilateral decoupling ---------------------------------------------------------

tau_hat <- matrix(1, nrow = J, ncol = J)
tau_hat[CHN_id, US_id] <- 1.01
tau_hat[US_id, CHN_id] <- 1.01

shock$gvc_chn_us <- shock$unchanged
shock$gvc_chn_us$tau_hat <- c(rep(as.vector(tau_hat), S * S), rep(1, J * J * S))

# decoupling us and eu from china ----------------------------------------------

tau_hat <- matrix(1, nrow = J, ncol = J)
tau_hat[CHN_id, US_id] <- 1.01
tau_hat[CHN_id, EU_id] <- 1.01
tau_hat[US_id, CHN_id] <- 1.01
tau_hat[EU_id, CHN_id] <- 1.01

shock$gvc_chn_us_with_eu <- shock$unchanged
shock$gvc_chn_us_with_eu$tau_hat <- c(rep(as.vector(tau_hat), S * S),
                                      rep(1, J * J * S))

# decoupling world from china gvc ----------------------------------------------

tau_hat <- matrix(1, nrow = J, ncol = J)
tau_hat[CHN_id, -CHN_id] <- 1.01
tau_hat[-CHN_id, CHN_id] <- 1.01

shock$gvc_chn_world <- shock$unchanged
shock$gvc_chn_world$tau_hat <- c(rep(as.vector(tau_hat), S * S),
                                 rep(1, J * J * S))

# decoupling world from china final --------------------------------------------

shock$final_chn_world <- shock$unchanged
shock$final_chn_world$tau_hat <- c(rep(1, J * J * S * S), rep(tau_hat, S))

# decoupling world from china total --------------------------------------------

shock$total_chn_world <- shock$unchanged
shock$total_chn_world$tau_hat <- rep(tau_hat, S * (S +1))

# one_way US-CHN ---------------------------------------------------------------

tau_hat <- matrix(1, nrow = J, ncol = J)
tau_hat[CHN_id, US_id] <- 1.01

shock$gvc_chn_us_oneway <- shock$unchanged
shock$gvc_chn_us_oneway$tau_hat <- c(rep(as.vector(tau_hat), S * S),
                                     rep(1, J * J * S))

# unilateral US ----------------------------------------------------------------

tau_hat <- matrix(1, nrow = J, ncol = J)
tau_hat[-US_id, US_id] <- 1.01

shock$gvc_us_unilateral <- shock$unchanged
shock$gvc_us_unilateral$tau_hat <- c(rep(as.vector(tau_hat), S * S),
                                     rep(1, J * J * S))

# no deficit shock -------------------------------------------------------------

shock$no_deficit <- within(shock$unchanged, varkappa_hat <- rep(0, J))

# Compute Counterfactual Baselines =============================================

# world GDP fixed throughout

# first run the no-deficit shock -----------------------------------------------

parameters <-
  list("epsilon" = trade_elasticities,
       "varphi" = intersec_elas,
       "mobility" = "perfect")

wiot <- mutate(
  wiot,
  no_deficit_cf = map(simulation_data, calc_cf, shock$no_deficit, parameters, nthreads = 8),
  # apply counterfactual changes
  no_deficit = map2(simulation_data, no_deficit_cf,
                    ~ within(.x, {
                      R <- R * .y$R_hat
                      pi <- pi * .y$pi_hat
                      D <- .y$D_prime
                    }))
)


data <-
  list(
    "wiot_2014_nozero" = filter(wiot, year == 2014, !includes_zeros) %>%
      pull(simulation_data) %>% .[[1]],
    "no_deficit" = filter(wiot, year == 2014,!includes_zeros) %>%
      pull(no_deficit) %>% .[[1]],
    "no_deficit_withzero" = filter(wiot, year == 2014, includes_zeros) %>%
      pull(no_deficit) %>% .[[1]],
    "no_deficit_2000" = filter(wiot, year == 2000,!includes_zeros) %>%
      pull(no_deficit) %>% .[[1]]
  )

shock$no_deficit <- NULL

# from these options calculate baselines for productivity shocks ---------------

# do all shocks for 2014 no_deficit data
baselines_cf <- tibble(data_name = "no_deficit",
                       shock_name = names(shock),
                       shock_size_percent = c(Inf),
                       mobility_type = "perfect")
# do gvc and gvc_except_eu shocks for varying strengths
baselines_cf <- tibble(data_name = "no_deficit",
                       shock_name = "gvc",
                       shock_size_percent = c(1, 5, 10, 25, 50, 100, 200),
                       mobility_type = "perfect") %>%
  bind_rows(baselines_cf)
baselines_cf <- tibble(data_name = "no_deficit",
                       shock_name = "gvc_except_eu",
                       shock_size_percent = c(1, 5, 10, 25, 50, 100, 200),
                       mobility_type = "perfect") %>%
  bind_rows(baselines_cf)

# calculate shocks from no_deficit baseline
baselines_cf <- baselines_cf %>%
  mutate(counterfactual =
    pmap(baselines_cf,
         function(data_name, shock_name, shock_size_percent, mobility_type) {
           print(paste0("Start data: ", data_name,
                        " shock: ", shock_name,
                        " mobility: ", mobility_type,
                        " shock size: ", shock_size_percent))
           calc_cf(data[[data_name]],
                   within(shock[[shock_name]],
                          tau_hat[tau_hat==1.01] <- 1 + shock_size_percent/100),
                   within(parameters, mobility <- mobility_type),
                   tol = 1e-7,
                   maxiter = 100000,
                   nthreads = 8)
           }))


# Set up baselines for future shocks
baselines <- baselines_cf %>%
  mutate(baseline = pmap(baselines_cf,
                         function(data_name, shock_name, shock_size_percent, mobility_type, counterfactual) {
                           baseline    <- data[[data_name]]
                           baseline$R  <- baseline$R  * counterfactual$R_hat
                           baseline$D  <- counterfactual$D_prime
                           baseline$pi <- baseline$pi * counterfactual$pi_hat
                           return(baseline)
                         }))

# Output to Stata --------------------------------------------------------------

baseline_to_df <- function(data_name, shock_name, shock_size_percent, mobility_type,
                           counterfactual, baseline) {
  S <- length(data[[data_name]]$sector_id)
  J <- length(data[[data_name]]$location_id)
  tibble(country = rep(data[[data_name]]$location_id, S),
         sector  = rep(data[[data_name]]$sector_id, each = J),
         Q_hat = as.vector(counterfactual$Q_hat),
         C_hat = rep(as.vector(counterfactual$C_hat), S),
         R_hat = as.vector(counterfactual$R_hat),
         P_rjU_hat = tail(counterfactual$P_hat, S*J),
         L_hat = as.vector(counterfactual$L_hat),
         L_tilde_hat = as.vector(counterfactual$L_tilde_hat),
         w_hat = R_hat / L_hat,
         alpha = as.vector(data[[data_name]]$alpha)) %>%
    group_by(country) %>%
    mutate(P_hat = prod(P_rjU_hat^alpha),
           realwage_hat = w_hat / P_hat)
}

baselines_stata <- baselines %>%
  mutate(counterfactual_df = pmap(baselines, baseline_to_df)) %>%
  select(-counterfactual, -baseline) %>%
  unnest(counterfactual_df)

# save and get rid of due to large memory usage
saveRDS(baselines, path$baselines)
rm(baselines)

# Calculate some baselines in the deficit_world and 0_trade_flows world, done
# separately due to memory/time constraints only

baselines_deficit_cf <- tibble(data_name = "wiot_2014_nozero",
                               shock_name = "gvc_except_eu",
                               shock_size_percent = c(10, Inf),
                               mobility_type = "perfect")
baselines_deficit_cf <- tibble(data_name = "no_deficit_withzero",
                               shock_name = "gvc_except_eu",
                               shock_size_percent = 10,
                               mobility_type = "perfect") %>%
  bind_rows(baselines_deficit_cf)
baselines_deficit_cf <- tibble(data_name = "wiot_2014_nozero",
                               shock_name = "gvc",
                               shock_size_percent = c(10, Inf),
                               mobility_type = "perfect") %>%
  bind_rows(baselines_deficit_cf)
baselines_deficit_cf <- tibble(data_name = "no_deficit_withzero",
                               shock_name = "gvc",
                               shock_size_percent = c(10, 200),
                               mobility_type = "perfect") %>%
  bind_rows(baselines_deficit_cf)


baselines_deficit_cf <- baselines_deficit_cf %>%
  mutate(counterfactual =
           pmap(baselines_deficit_cf,
                function(data_name, shock_name, shock_size_percent,
                         mobility_type) {
                  if(data_name == "wiot_2014_nozero" & is.infinite(shock_size_percent)) {
                    method = "direct"
                  } else {
                    method = "automatic"
                  }
                  print(paste0("Start data: ", data_name,
                               " shock: ", shock_name,
                               " mobility: ", mobility_type,
                               " shock size: ", shock_size_percent))
                  calc_cf(data[[data_name]],
                          within(shock[[shock_name]],
                                 tau_hat[tau_hat == 1.01] <-
                                   1 + shock_size_percent / 100),
                          within(parameters,
                                 mobility <- mobility_type),
                          tol = 1e-7,
                          maxiter = 100000,
                          method = method,
                          nthreads = 8)}))

# Set up baselines for future shocks
baselines_deficit <- baselines_deficit_cf %>%
  mutate(baseline = pmap(baselines_deficit_cf,
                         function(data_name, shock_name, shock_size_percent,
                                  mobility_type, counterfactual) {
                           baseline    <- data[[data_name]]
                           baseline$R  <- baseline$R  * counterfactual$R_hat
                           baseline$D  <- counterfactual$D_prime
                           baseline$pi <- baseline$pi * counterfactual$pi_hat
                           return(baseline)
                         }))

# Output to Stata

baselines_stata <- baselines_deficit %>%
  mutate(counterfactual_df = pmap(baselines_deficit, baseline_to_df)) %>%
  select(-counterfactual, -baseline) %>%
  unnest(counterfactual_df) %>%
  bind_rows(baselines_stata)

saveRDS(baselines_deficit, path$baselines_deficit)
rm(baselines_deficit)

# Policy Robustness Checks -----------------------------------------------------

# To redo table 1 for different shock sizes
baselines_cf <- tibble(data_name = "no_deficit",
                       shock_name = "gvc_chn_us",
                       shock_size_percent = c(10, 50, 100, 200),
                       mobility_type = "perfect")
baselines_cf <- tibble(data_name = "no_deficit",
                       shock_name = "gvc_chn_us_oneway",
                       shock_size_percent = c(10, 50, 100, 200),
                       mobility_type = "perfect") %>%
  bind_rows(baselines_cf)
baselines_cf <- tibble(data_name = "no_deficit",
                       shock_name = "gvc_us_unilateral",
                       shock_size_percent = c(10, 50, 100, 200),
                       mobility_type = "perfect") %>%
  bind_rows(baselines_cf)
baselines_cf <- tibble(data_name = "no_deficit",
                       shock_name = "gvc_chn_us_with_eu",
                       shock_size_percent = c(10, 50, 100, 200),
                       mobility_type = "perfect") %>%
  bind_rows(baselines_cf)

# calculate shocks from no_deficit baseline
baselines_cf <- baselines_cf %>%
  mutate(counterfactual =
           pmap(baselines_cf,
                function(data_name, shock_name, shock_size_percent,
                         mobility_type) {
                  print(paste0("Start data: ", data_name,
                               " shock: ", shock_name,
                               " mobility: ", mobility_type,
                               " shock size: ", shock_size_percent))
                  calc_cf(data[[data_name]],
                          within(shock[[shock_name]],
                                 tau_hat[tau_hat == 1.01] <-
                                   1 + shock_size_percent / 100),
                          within(parameters,
                                 mobility <- mobility_type),
                          tol = 1e-7,
                          maxiter = 100000,
                          nthreads = 8)}))


# Set up baselines for future shocks
baselines_policy <- baselines_cf %>%
  mutate(baseline = pmap(baselines_cf,
                         function(data_name, shock_name, shock_size_percent, mobility_type, counterfactual) {
                           baseline    <- data[[data_name]]
                           baseline$R  <- baseline$R  * counterfactual$R_hat
                           baseline$D  <- counterfactual$D_prime
                           baseline$pi <- baseline$pi * counterfactual$pi_hat
                           return(baseline)
                         }))

baselines_stata <- baselines_policy %>%
  mutate(counterfactual_df = pmap(baselines_policy, baseline_to_df)) %>%
  select(-counterfactual, -baseline) %>%
  unnest(counterfactual_df) %>%
  bind_rows(baselines_stata)

saveRDS(baselines_policy, path$baselines_policy)
rm(baselines_policy)

# Save =========================================================================

write_dta(baselines_stata, path$baselines_stata)