In VKudlay/GEGravity: Solver for Armington-CES General Equilibrium Trade Model

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(GEGravity)
if (!requireNamespace("alpaca", quietly = TRUE)) {
  stop("Package \"alpaca\" needed for about.Rmd and example.Rmd to work. Please install it.",
    call. = FALSE)
}
if (!requireNamespace("data.table", quietly = TRUE)) {
  stop("Package \"data.table\" needed for about.Rmd and example.Rmd to work. Please install it.",
    call. = FALSE)
}

This example RMD uses the example data TradeData0014, which is included with this package. The data set consists of a panel of 44 countries trading with one another over the years 2000-2014. The trade data uses aggregated trade flows based on WIOD and information on FTAs is from the NSF-Kellogg database maintained by Scott Baier and Jeff Bergstrand. To find out more about the trade data, see "Timmer, Dietzenbacher, Los, Stehrer, and de Vries (2015).

data(TradeData0014) # loads data included with package
head(TradeData0014)

Suppose the researcher wishes to use this data set to quantify general equilibrium trade and welfare effects of the EU enlargements that took place between 2000-2014. To first obtain the partial effects of these enlargements on trade flows, a PPML (Poisson pseudo-maximum likelihood) gravity specification may be used. Specifically, to obtain the "partial" estimates of the effects of EU enlargements on trade, we can use the following three-way gravity specification:

$$X_{ijt} = \exp\bigg[\ln(\alpha_{it}) + \ln(\alpha_{jt}) + \ln(\alpha_{ij}) + \beta \cdot FTA_{ijt}\bigg] + e_{ijt}$$

The Stata equivalent of this package uses the ppmlhdfe command created by Correia, Guimarães, & Zylkin (2019) to achieve this. In our example, we will use the 'alpaca' library by Amrei Stammann to do the same via the \code{feglm} (Fixed-Effect Generalized Linear Model) function.

library(alpaca)  # Needed for partial coefficients

# Generate foreign trade subset
f_trade <- TradeData0014[TradeData0014$exporter != TradeData0014$importer,]

# classify FEs for components to be absorbed (finding variable interactions)
f_trade$exp_year <- interaction(f_trade$expcode, f_trade$year)
f_trade$imp_year <- interaction(f_trade$impcode, f_trade$year)
f_trade$pair     <- interaction(f_trade$impcode, f_trade$expcode)

# Fit gravity model using PPML with exporter-time, importer-time, and exporter-importer FEs
partials <- feglm(
  formula = trade ~ eu_enlargement + other_fta | exp_year + imp_year + pair,
  data    = f_trade,
  family  = poisson()
)$coefficient  # We just need the coefficients for computation

print(round(partials, 3))

In addition to estimating the effects of EU enlargements on new EU pairs, this example also controls for any other FTAs signed during the period. Each of these variables is coded as a dummy variable that becomes 1 when the agreement goes into effect for a given pair. The estimated coefficient for \code{eu_enlargements} is 0.224, implying that the expansion of the EU had an average partial effect of $e^{0.224} − 1 = 25.1\%$ on trade between new EU members and existing members. With clustered standard errors, this estimate is statistically significant at the $p < .01$ significance level.

Before proceeding, some further pre-processing is needed. Specifically, we need to supply the function with the partial effects of joining the EU estimated above for each of the appropriate pairs. To do this, we create a new variable whose value is either the eu_enlargement partial computed above or 0. We will use the year 2000 as the baseline for the general equilibrium counterfactual, so we assign these partial effects in the year 2000.

For Stata users, the Stata equivalent of this next step is:

sort exporter importer year
by exporter importer: gen new_eu_pair = (eu_enlargement[_N]-eu_enlargement[1])                   
by exporter importer: gen eu_effect = _b[eu_enlargement] * new_eu_pair

In R, we do the following:

# Sort matrix to make it easier to find imp/exp pairs
t_trade <- TradeData0014[order(
  TradeData0014$exporter,
  TradeData0014$importer,
  TradeData0014$year),
]

t_trade$new_eu_pair <- NA
t_trade$eu_effect   <- NA   # this creates a new column that will contain partial effect of EU membership for new EU pairs
i <- 1
# Effect of EU entrance on country based on partial, if entry happened
invisible(by(t_trade, list(t_trade$expcode, t_trade$impcode), function(row) {
  # Was a new EU pair created within time span?
  t_trade[i:(i+nrow(row)-1), "new_eu_pair"] <<- diff(row$eu_enlargement, lag=nrow(row)-1)
  i <<- i + nrow(row)
}))

# If added to EU, give it the computed partial eu_enlargement coefficient (0.224) as the effect
t_trade$eu_effect = t_trade$new_eu_pair * partials[1]

To finalize the data for the counterfactual, we will use the year 2000 as the baseline year. The data we will feed to the 'ge_gravity' command looks like this:

# Data to be finally fed to the function (we base the counterfactual on the year 2000.)
ge_baseline_data <- t_trade[t_trade$year == 2000,]   # In example, 1892 Entries, 5676 removed

# head(data) # First 10 rows
ge_baseline_data[sample(1:nrow(ge_baseline_data), 10),]  # 10 random rows

Now that we have the data we need, we can run the 'ge_gravity' function as follows:

ge_results <- ge_gravity(
  exp_id = ge_baseline_data$expcode,    # Origin country associated with each observation
  imp_id = ge_baseline_data$impcode,    # Destination country associated with each observation
  flows  = ge_baseline_data$trade,      # Observed trade flows in the data for the year being used as the baseline
  beta   = ge_baseline_data$eu_effect,  # “Partial” change in trade, obtained as coefficient from gravity estimation
  theta  = 4,               # Trade elasticity
  mult   = FALSE,           # Assume trade balance is an additive component of national expenditure
  data   = ge_baseline_data
)
ge_results[sample(1:nrow(ge_results), 10),c(1:2,5:6,12:16)] # 10 random rows

This assumes a standard trade elasticity value of $\theta = 4$. The input for $\beta$ is given by the variable we created called eu_effect, which is equal to 0.224 for new EU pairs formed during the period and equal to 0 otherwise. Because of the small size of the sample, it solves almost instantly. Sample results for the first 10 rows in the data are shown above for the counterfactual trade level and the associated changes in welfare, real wages, nomimal wages, and the local price index. Click on the right arrow to scroll to the right if they are not all displayed above. For the latter four variables, the number shown is the result computed for the exporting country. Unsurprisingly, the new EU members (Bulgaria, Croatia, Czech Republic, Estonia, Hungary, Latvia, Lithuania, Malta, Poland, Romania, Slovakia, and Slovenia) realize the largest welfare gains from their joining the EU, with existing EU countries also gaining. All countries not included in the agreement experience small losses due to trade diversion, with the largest losses accruing to Russia.

We can also change how trade imbalances enter the model. The default (exibited above with mult = FALSE) is to assume that they enter expenditure additively (i.e., $E_j = Y_j + D_j$), but one can also change the model so that expenditure is instead a fixed multiple of income (i.e., let $E_j = \delta_j Y_j$.) by setting mult = TRUE. While using multiplicative imbalances instead of additive balances changes the results slightly, they are still qualitatively very similar.

An important point about the above exercises is that the initial partial effect is estimated with some uncertainty. The GE results that were calculated may paint a misleading picture because they do not take this uncertainty into account. For this reason, it is considered good practice to use a bootstrap method to construct confidence intervals for the GE calculations. This type of procedure can be made easier using ge_gravity function:

library(data.table)  # needed below for bootstrap procedure

# Helper function for shuffling the data with replacement.
# This allows us to shuffle by pair rather than treating each observation as independent.
get_bootdata <- function(data=list(),id) {
  uniq_id <- levels(factor(data[,id]))
  draw    <- sort(sample(uniq_id,replace=TRUE))
  draw  <- data.table(draw)[,.N,by=draw]  # count duplicate pairs in bootstrap sample
  colnames(draw)[1] <- id
  boot_data  <- merge(data,draw,"pair",all.x=FALSE,all.y=FALSE)
  boot_index <- rep(row.names(boot_data), boot_data$N)  # replicates pairs drawn multiple times after merge
  boot_data  <- boot_data[matrix(boot_index), ]
  boot_data$rep <- (as.numeric(rownames(boot_data)) %% 1)*10+1
  return(data.frame(boot_data))
}

# set up for pair bootstrap
set.seed(12345)
bootreps      <- 20
TradeData0014[,"pair"] <- interaction(TradeData0014$expcode,TradeData0014$impcode) # This is the ID we will use for resampling.

# Initialize matrices for saving results
x             <- TradeData0014[,c("eu_enlargement","other_fta")]
save_partials <- matrix(nrow = bootreps, ncol = ncol(x), dimnames = list(1:bootreps,colnames(x)))
GE_effects    <- ge_results[,11:15]
save_GE       <- matrix(nrow = bootreps, ncol = ncol(GE_effects), dimnames = list(1:bootreps,colnames(GE_effects)))

# generate bootstrapped gravity estimates using alpaca's feglm function (20 boot reps)
library(alpaca)
for (b in 1:bootreps) {

    # This step shuffles the data using the get_bootdata() function defined above
    boot_data <- get_bootdata(TradeData0014,"pair")

    # These next few steps are exactly the same as the ones we used above to estimate the partial effects

    # Generate foreign trade subset
    f_trade <- boot_data[boot_data$exporter != boot_data$importer,]

    # Normalize trade data to unit interval
    f_trade$trade <- f_trade$trade / max(f_trade$trade)

    # classify FEs to be absorbed
    f_trade$exp_year <- interaction(f_trade$expcode, f_trade$year)
    f_trade$imp_year <- interaction(f_trade$impcode, f_trade$year)
    f_trade$pair     <- interaction(f_trade$impcode, f_trade$expcode)

    # Estimate and save partial effects
    save_partials[b,] <- feglm(
    formula = trade ~ eu_enlargement + other_fta | exp_year + imp_year + pair,
    data    = f_trade,
    family  = poisson()
    )$coefficient  # We just need the coefficients for computation
}


# Obtain bootstrapped GE results based on bootstrapped partial effects
bootstrap_GE_results <- ge_baseline_data
for (b in 1:bootreps) {

  # set up baseline data using the estimated partial effect from bootstrap b
  boot_ge_data <- ge_baseline_data
  boot_ge_data$eu_effect <- save_partials[b,1] * boot_ge_data$new_eu_pair 

  # run GE_gravity
  temp <- ge_gravity(
  exp_id = boot_ge_data$expcode,    # Origin country associated with each observation
  imp_id = boot_ge_data$impcode,    # Destination country associated with each observation
  flows  = boot_ge_data$trade,      # Observed trade flows in the data for the year being used as the baseline
  beta   = boot_ge_data$eu_effect,  # “Partial” change in trade, obtained as coefficient from gravity estimation
  theta  = 4,               # Trade elasticity
  mult   = FALSE,           # Assume trade balance is an additive component of national expenditure
  data   = boot_ge_data
  )

  # store results
  bootstrap_GE_results[,paste0("welf",b)] <- temp[,"welfare"]
  bootstrap_GE_results[,paste0("trade",b)] <- temp[,"new_trade"]
}

# get bootstrapped means, SDs, and 95% CIs for partial effects
colMeans(save_partials)
apply(save_partials, 2, sd)
apply(save_partials, 2, function(x) quantile(x, probs = .975))
apply(save_partials, 2, function(x) quantile(x, probs = .025))

# Get 95% CIs for GE effects
temp <- bootstrap_GE_results[,12:51]
temp <- temp[,order(colnames(temp))]

bootstrap_GE_results[,"lb_welf"] <- apply(temp[,21:40], 1, function(x) quantile(x, probs = .025))
bootstrap_GE_results[,"ub_welf"] <- apply(temp[,21:40], 1, function(x) quantile(x, probs = .975))
bootstrap_GE_results[,"lb_trade"] <- apply(temp[,1:20], 1, function(x) quantile(x, probs = .025))
bootstrap_GE_results[,"ub_trade"] <- apply(temp[,1:20], 1, function(x) quantile(x, probs = .975))

disp_cols <- c("exporter","importer","year","trade","lb_welf","ub_welf","lb_trade","ub_trade")
bootstrap_GE_results[sample(1:nrow(ge_results), 10),disp_cols] # 10 random rows; GE welfare effects are for *exporter*