Chapter 3 General equilibrium trade policy analysis with structural gravity

3.1 Trade without borders

3.1.1 Initial data

Unlike the previous chapter, we shall proceed by alternating both data transforming and regressions. In the previous chapter it was possible to first process the datasets and then fit the regressions, but here we need the regressions’ output in order to create new variables. In any case we will follow quite similar steps to the last chapter.

To do what is shown in box #1 from page 104 in Yotov et al. (2016), we need to convert “DEU” in both exporter and importer columns to “0-DEU.” We could have used “AAA,” the book uses “ZZZ” but in R “ZZZ” won’t be treated as the reference factor. It is important to mention that box #1 doesn’t show a previous step that is mentioned in page 103, which is to filter and keep observations for the year 2006 only.

ch2_application1 <- yotov_data("ch2_application1") %>% 
  filter(year == 2006) %>% 
  mutate(
    log_dist = log(dist),
    intl = ifelse(exporter != importer, 1, 0),
    exporter = ifelse(exporter == "DEU", "0-DEU", exporter),
    importer = ifelse(importer == "DEU", "0-DEU", importer)
  ) %>% 
  
  # Create Yit
  group_by(exporter, year) %>% 
  mutate(y = sum(trade)) %>% 
  
  # Create Eit
  group_by(importer, year) %>% 
  mutate(e = sum(trade)) %>% 
  
  # Create Er
  ungroup() %>% 
  mutate(e_r = max(ifelse(importer == "0-DEU", e, NA), na.rm = T))

What we did at this point is exactly the same as the previous chapter, we are summarising our data and adding new columns in terms of descriptive statistics (i.e. the maximum for e_r), and this is our dataset up to this point:

ch2_application1
## # A tibble: 4,761 x 14
##    exporter importer pair_id  year   trade   dist  cntg  lang  clny log_dist
##    <chr>    <chr>      <int> <int>   <dbl>  <dbl> <int> <int> <int>    <dbl>
##  1 ARG      ARG        12339  2006 32313.    534.     0     0     0     6.28
##  2 ARG      AUS            1  2006   108.  12045.     0     0     0     9.40
##  3 ARG      AUT            2  2006    25.2 11751.     0     0     0     9.37
##  4 ARG      BEL            4  2006   189.  11305.     0     0     0     9.33
##  5 ARG      BGR            3  2006    17.0 12116.     0     0     0     9.40
##  6 ARG      BOL            6  2006   392.   1866.     1     1     0     7.53
##  7 ARG      BRA            8  2006  6198.   2392.     1     0     0     7.78
##  8 ARG      CAN           10  2006   339.   9391.     0     0     0     9.15
##  9 ARG      CHE           12  2006   278.  11233.     0     0     0     9.33
## 10 ARG      CHL           14  2006  2456.   1157.     1     1     0     7.05
## # … with 4,751 more rows, and 4 more variables: intl <dbl>, y <dbl>, e <dbl>,
## #   e_r <dbl>

3.1.2 Step I: Solve the baseline model

We start by fitting the next model:

\[ \begin{align} X_{ij,t} =& \:\exp\left[\pi_{i,t} + \chi_{i,t} + \beta_1 \log(DIST)_{i,j} + \beta_2 CNTG_{i,j} + \beta_3 INTL_{i,j}\right] \times \varepsilon_{ij,t} \end{align} \]

With the data from above, the model specification is straightforward:

fit_baseline_app1 <- glm(
  trade ~ 0 + log_dist + cntg + intl + exporter + importer,
  family = quasipoisson(link = "log"),
  data = ch2_application1
)

For now, we will concentrate on the fitted values and shall ignore the clustered standard errors in the next paragraphs, but still we can show the robust estimation by using the yotov_clustered_glm() function:

yotov_clustered_glm(fit_baseline_app1$formula, ch2_application1)
## 
## z test of coefficients:
## 
##           Estimate Std. Error  z value  Pr(>|z|)    
## log_dist -0.791288   0.061857 -12.7922 < 2.2e-16 ***
## cntg      0.673646   0.137746   4.8905 1.006e-06 ***
## intl     -2.474450   0.147070 -16.8250 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

With the estimated model, we can proceed as in box #1 from page 105 in Yotov et al. (2016) in order to construct the variables for export and import fixed effects:

ch2_application1 <- ch2_application1 %>% 
  left_join(
    yotov_fixed_effects(fit_baseline_app1),
    c("exporter", "importer")
  )

Still following box #1, we need to compute the variables of bilateral trade costs and multilateral resistances:

ch2_application1 <- ch2_application1 %>% 
  mutate(
    tij_bln = exp(fit_baseline_app1$coefficients["log_dist"] * log_dist +
                  fit_baseline_app1$coefficients["cntg"] * cntg +
                  fit_baseline_app1$coefficients["intl"] * intl),
  
    # outward multilateral resistance (omr)
    omr_bln = y * (e_r / exp(fe_exporter)),
    
    # inward multilateral resistance (imr)
    imr_bln = e / (exp(fe_importer) * e_r)
  )

To complete this stage of the estimation, we need to create a column with the estimated international trade for given output and expenditures. We start by adding a column, tradehat_bln, with the regression output, and then we group by exporter and summarise to obtain the required column xi_bln:

ch2_application1 <- ch2_application1 %>% 
  mutate(tradehat_bln = predict(fit_baseline_app1, ch2_application1, "response")) %>% 
  group_by(exporter) %>% 
  mutate(xi_bln = sum(tradehat_bln * (exporter != importer))) %>% 
  ungroup()

The dataset at this point is not much different to the previous chapter, with the difference that xi_bln depends on the regression output, so now we are using statistical inference to create new variables which depend on the existing ones. We can give a look to these new variables and the ones that create those:

ch2_application1 %>% 
  select(xi_bln, tradehat_bln, trade, log_dist, cntg, intl, exporter, importer)
## # A tibble: 4,761 x 8
##    xi_bln tradehat_bln   trade log_dist  cntg  intl exporter importer
##     <dbl>        <dbl>   <dbl>    <dbl> <int> <dbl> <chr>    <chr>   
##  1 34709.      24853.  32313.      6.28     0     0 ARG      ARG     
##  2 34709.        862.    108.      9.40     0     1 ARG      AUS     
##  3 34709.        172.     25.2     9.37     0     1 ARG      AUT     
##  4 34709.        344.    189.      9.33     0     1 ARG      BEL     
##  5 34709.         55.5    17.0     9.40     0     1 ARG      BGR     
##  6 34709.        222.    392.      7.53     1     1 ARG      BOL     
##  7 34709.       6355.   6198.      7.78     1     1 ARG      BRA     
##  8 34709.       1121.    339.      9.15     0     1 ARG      CAN     
##  9 34709.        260.    278.      9.33     0     1 ARG      CHE     
## 10 34709.       3084.   2456.      7.05     1     1 ARG      CHL     
## # … with 4,751 more rows

3.1.3 Step II: Define a counterfactual scenario

Box #2 from page 105 in Yotov et al. (2016) proposes two alternatives to define counterfactual scenario of removing international borders.

The first alternative is to eliminate the border variable and then generate the logged trade costs used in the constraint:

ch2_application1 <- ch2_application1 %>% 
  mutate(
    tij_cfl = exp(fit_baseline_app1$coefficients["log_dist"] * log_dist +
                  fit_baseline_app1$coefficients["cntg"] * cntg),
    log_tij_cfl = log(tij_cfl)
  )

The second alternative (which won’t be used here, but is left as an exercise) is to define a new counterfactual border variable:

ch2_application1 <- ch2_application1 %>% 
  mutate(
    intl_cfl = 0,
    tij_cfl = exp(fit_baseline_app1$coefficients["log_dist"] * log_dist +
                  fit_baseline_app1$coefficients["cntg"] * cntg +
                  fit_baseline_app1$coefficients["intl"] * intl_cfl),
    log_tij_cfl = log(tij_cfl)
  )

Again, it’s convenient to explore the column that we just created:

ch2_application1 %>% 
  select(tij_cfl, log_dist, cntg)
## # A tibble: 4,761 x 3
##     tij_cfl log_dist  cntg
##       <dbl>    <dbl> <int>
##  1 0.00695      6.28     0
##  2 0.000590     9.40     0
##  3 0.000602     9.37     0
##  4 0.000620     9.33     0
##  5 0.000587     9.40     0
##  6 0.00506      7.53     1
##  7 0.00416      7.78     1
##  8 0.000718     9.15     0
##  9 0.000624     9.33     0
## 10 0.00739      7.05     1
## # … with 4,751 more rows

3.1.4 Step III: Solve the counterfactual model

We need to fit a model similar to the model from step I, the constrained gravity model, where \(\pi_{j,t}\) and \(\chi_{j,t}\) are altered:

\[ \begin{align} X_{ij,t} =& \:\exp\left[\pi_{i,t}^{CFL} + \chi_{i,t}^{CFL} + \beta_1 \log(DIST)_{i,j} + \beta_2 CNTG_{i,j} + \beta_3 INTL_{i,j}\right] \times \varepsilon_{ij,t} \end{align} \]

Box #1 from page 106 in Yotov et al. (2016) estimates the constrained gravity model with the PPML estimator by using an offset argument, and this is straightforward in R:

fit_counterfactual_app1 <- glm(
  trade ~ 0 + exporter + importer + offset(log_tij_cfl),
  family = quasipoisson(link = "log"),
  data = ch2_application1
)

Unlike step I, to construct the variables for export and import fixed effects, we’ll obtain the variables fe_exporter.x and fe_exporter.y because we already added an fe_exporter column to the dataset after obtaining the baseline model.

dplyr is wise enough to rename the two fe_exporter columns (the same applies to the imports) but we need to specify that we want to join by exporter and importer instead of all columns with shared names. In any case, it is better to rename those columns and provide informative names:

ch2_application1 <- ch2_application1 %>% 
  left_join(
    yotov_fixed_effects(fit_counterfactual_app1),
    by = c("exporter", "importer")
  ) %>% 
  rename(
    fe_exporter_bln = fe_exporter.x,
    fe_exporter_cfl = fe_exporter.y,
    fe_importer_bln = fe_importer.x,
    fe_importer_cfl = fe_importer.y
  )

Now we go for Box #2 from page 106 in Yotov et al. (2016) where the variables of bilateral trade costs and multilateral resistances are obtained:

ch2_application1 <- ch2_application1 %>% 
  mutate(
    # outward multilateral resistance (omr)
    omr_cfl = y * (e_r / exp(fe_exporter_cfl)),

    # inward multilateral resistance (imr)
    imr_cfl = e / (exp(fe_importer_cfl) * e_r)
  )

Box #2 also shows how to compute the conditional general equilibrium effects of trade. This is very similar to what we did in step I:

ch2_application1 <- ch2_application1 %>% 
  mutate(tradehat_cfl = predict(fit_counterfactual_app1, ch2_application1, "response")) %>% 
  group_by(exporter) %>% 
  mutate(xi_cfl = sum(tradehat_cfl * (exporter != importer))) %>% 
  ungroup()

Box #1 from page 107 in Yotov et al. (2016) can be considerably simplified with R code. To construct the iterative procedure to converge to full endowment general equilibrium effects, we start by creating the required columns and parameters, so we will deviate from the original approach.

We start computing change in bilateral trade costs (change_tij) and trade deficit or surplus (phi):

# set the criteria of convergence
# taken from the literature (see the Stata code)
sigma <- 7

ch2_application1 <- ch2_application1 %>% 
  mutate(
    change_tij = tij_cfl / tij_bln,
    phi = ifelse(importer == exporter, e / y, 0)
  ) %>% 
  group_by(exporter) %>% 
  mutate(phi = max(phi)) %>% 
  ungroup()

Now we compute change in prices for exporters (change_p_i) and importers (change_p_j):

ch2_application1 <- ch2_application1 %>% 
  group_by(exporter) %>% 
  mutate(change_p_i = ((exp(fe_exporter_cfl) / e_r) / (exp(fe_exporter_bln) / e_r))^(1 /(1 - sigma))) %>% 
  ungroup() %>% 
  
  group_by(importer) %>% 
  mutate(
    change_p_j = ifelse(importer == exporter, change_p_i, 0),
    change_p_j = max(change_p_j)
  ) %>% 
  ungroup()

Next, we need to compute the counterfactual trade flows (trade_cfl):

ch2_application1 <- ch2_application1 %>% 
  mutate(trade_cfl = tradehat_cfl * change_p_i * change_p_j)

Before going further, we have obtained the next dataset:

ch2_application1
## # A tibble: 4,761 x 34
##    exporter importer pair_id  year   trade   dist  cntg  lang  clny log_dist
##    <chr>    <chr>      <int> <int>   <dbl>  <dbl> <int> <int> <int>    <dbl>
##  1 ARG      ARG        12339  2006 32313.    534.     0     0     0     6.28
##  2 ARG      AUS            1  2006   108.  12045.     0     0     0     9.40
##  3 ARG      AUT            2  2006    25.2 11751.     0     0     0     9.37
##  4 ARG      BEL            4  2006   189.  11305.     0     0     0     9.33
##  5 ARG      BGR            3  2006    17.0 12116.     0     0     0     9.40
##  6 ARG      BOL            6  2006   392.   1866.     1     1     0     7.53
##  7 ARG      BRA            8  2006  6198.   2392.     1     0     0     7.78
##  8 ARG      CAN           10  2006   339.   9391.     0     0     0     9.15
##  9 ARG      CHE           12  2006   278.  11233.     0     0     0     9.33
## 10 ARG      CHL           14  2006  2456.   1157.     1     1     0     7.05
## # … with 4,751 more rows, and 24 more variables: intl <dbl>, y <dbl>, e <dbl>,
## #   e_r <dbl>, fe_exporter_bln <dbl>, fe_importer_bln <dbl>, tij_bln <dbl>,
## #   omr_bln <dbl>, imr_bln <dbl>, tradehat_bln <dbl>, xi_bln <dbl>,
## #   tij_cfl <dbl>, log_tij_cfl <dbl>, fe_exporter_cfl <dbl>,
## #   fe_importer_cfl <dbl>, omr_cfl <dbl>, imr_cfl <dbl>, tradehat_cfl <dbl>,
## #   xi_cfl <dbl>, change_tij <dbl>, phi <dbl>, change_p_i <dbl>,
## #   change_p_j <dbl>, trade_cfl <dbl>

To conclude the steps from Box #1 we need a while() loop and iterate until a convergence is reached. We need to duplicate some columns under new names for the loop operations, because we will overwrite them when using the iterative steps:

ch2_application1 <- ch2_application1 %>% 
  mutate(
    omr_cfl_0 = omr_cfl,
    imr_cfl_0 = imr_cfl,
    change_imr_full_0 = 1,
    change_omr_full_0 = 1,
    change_p_i_0 = change_p_i,
    change_p_j_0 = change_p_j,
    fe_exporter_cfl_0 = fe_exporter_cfl,
    fe_importer_cfl_0 = fe_importer_cfl,
    tradehat_0 = tradehat_cfl,
    e_r_cfl_0 = e_r
  )

Observe that we have duplicated columns now, let’s explore a few columns:

ch2_application1 %>% 
  select(omr_cfl, omr_cfl_0, change_imr_full_0, change_omr_full_0)
## # A tibble: 4,761 x 4
##    omr_cfl omr_cfl_0 change_imr_full_0 change_omr_full_0
##      <dbl>     <dbl>             <dbl>             <dbl>
##  1  36795.    36795.                 1                 1
##  2  36795.    36795.                 1                 1
##  3  36795.    36795.                 1                 1
##  4  36795.    36795.                 1                 1
##  5  36795.    36795.                 1                 1
##  6  36795.    36795.                 1                 1
##  7  36795.    36795.                 1                 1
##  8  36795.    36795.                 1                 1
##  9  36795.    36795.                 1                 1
## 10  36795.    36795.                 1                 1
## # … with 4,751 more rows

Here omr_cfl_0 will the consecutively updated with the forthcoming iterative process, and therefore on later iterations it will be different from it’s initial values we are seeing here. The change_ columns will also change.

Now we run the loop, which cannot be divided into smaller pieces because the \(N^{th}\) step depends on the \(1,2,\ldots,N-1\) previous steps:

# set parameters
max_dif <- 1
sd_dif <- 1
change_price_i_old <- 0

i <- 1
while(sd_dif > 1e-5 | max_dif > 1e-5) {
  ch2_application1 <- ch2_application1 %>% 
    mutate(trade_1 = tradehat_0 * change_p_i_0 * change_p_j_0 / (change_omr_full_0 * change_imr_full_0))
  
  # repeat the counterfactual model
  fit_counterfactual_app1_2 <- glm(
    trade_1 ~ 0 + exporter + importer + offset(log_tij_cfl),
    family = quasipoisson(link = "log"),
    data = ch2_application1
  )
  
  ch2_application1 <- ch2_application1 %>% 
    left_join(
      yotov_fixed_effects(fit_counterfactual_app1_2),
      by = c("exporter", "importer")
    )
  
  # compute the conditional general equilibrium effects of trade
  ch2_application1 <- ch2_application1 %>% 
    mutate(tradehat_1 = predict(fit_counterfactual_app1_2, ch2_application1, "response")) %>% 
    group_by(exporter) %>% 
    mutate(y_cfl_1 = sum(tradehat_1)) %>% 
    ungroup() %>% 
    
    mutate(e_cfl_1 = ifelse(importer == exporter, phi * y_cfl_1, 0)) %>% 
    group_by(importer) %>% 
    mutate(e_cfl_1 = max(e_cfl_1)) %>% 
    ungroup() %>% 
    
    mutate(
      e_r_cfl_1 = ifelse(importer == "0-DEU", e_cfl_1, 0),
      e_r_cfl_1 = max(e_r_cfl_1)
    )
  
  # compute the change in prices for exporters and importers
  ch2_application1 <- ch2_application1 %>% 
    mutate(change_p_i_1 = ((exp(fe_exporter) / e_r_cfl_1) / 
      (exp(fe_exporter_cfl_0) / e_r_cfl_0))^(1 / (1 - sigma)))
  
  # compute the change in prices for exporters and importers
  ch2_application1 <- ch2_application1 %>% 
    group_by(importer) %>% 
    mutate(
      change_p_j_1 = ifelse(importer == exporter, change_p_i_1, 0),
      change_p_j_1 = max(change_p_j_1)
    ) %>% 
    ungroup()
  
  # compute both outward and inward multilateral resistance
  ch2_application1 <- ch2_application1 %>% 
    mutate(
      omr_cfl_1 = (y_cfl_1 * e_r_cfl_1) / exp(fe_exporter),
      imr_cfl_1 = e_cfl_1 / (exp(fe_importer) * e_r_cfl_1)
    )
  
  # update the differences
  max_dif <- abs(max(ch2_application1$change_p_i_0 - change_price_i_old))
  sd_dif <- sd(ch2_application1$change_p_i_0 - change_price_i_old)
  change_price_i_old <- ch2_application1$change_p_i_0
  
  # compute changes in outward and inward multilateral resistance
  ch2_application1 <- ch2_application1 %>% 
    mutate(
      change_omr_full_1 = omr_cfl_1 / omr_cfl_0,
      change_imr_full_1 = imr_cfl_1 / imr_cfl_0,
      omr_cfl_0 = omr_cfl_1,
      imr_cfl_0 = imr_cfl_1,
      change_omr_full_0 = change_omr_full_1,
      change_imr_full_0 = change_imr_full_1,
      change_p_i_0 = change_p_i_1,
      change_p_j_0 = change_p_j_1,
      fe_exporter_cfl_0 = fe_exporter,
      fe_importer_cfl_0 = fe_importer,
      tradehat_0 = tradehat_1,
      e_r_cfl_0 = e_r_cfl_1
    ) %>% 
    select(-fe_exporter, -fe_importer)
  
  i <- i + 1
}

As an example, see how some columns changes:

ch2_application1 %>% 
  select(omr_cfl, omr_cfl_0, change_imr_full_0, change_omr_full_0)
## # A tibble: 4,761 x 4
##    omr_cfl omr_cfl_0 change_imr_full_0 change_omr_full_0
##      <dbl>     <dbl>             <dbl>             <dbl>
##  1  36795.    54561.              1.00              1.43
##  2  36795.    54561.              1.00              1.43
##  3  36795.    54561.              1.00              1.43
##  4  36795.    54561.              1.00              1.43
##  5  36795.    54561.              1.00              1.43
##  6  36795.    54561.              1.00              1.43
##  7  36795.    54561.              1.00              1.43
##  8  36795.    54561.              1.00              1.43
##  9  36795.    54561.              1.00              1.43
## 10  36795.    54561.              1.00              1.43
## # … with 4,751 more rows

Box #1 from page 108 in Yotov et al. (2016) shows the steps to obtain different endowments, which can be divided into smaller pieces.

We start computing the full endowment general equilibrium of factory-gate price (change_p_i_full and change_p_j_full) and the full endowment general equilibrium of output (y_full):

ch2_application1 <- ch2_application1 %>% 
  mutate(
    change_p_i_full = ((exp(fe_exporter_cfl_0) / e_r_cfl_0) /
                         (exp(fe_exporter_bln) / e_r))^(1 / (1 - sigma)),
    change_p_j_full = change_p_i_full * (exporter == importer)
  ) %>% 
  group_by(importer) %>% 
  mutate(change_p_j_full = max(change_p_j_full)) %>% 
  ungroup() %>% 
  mutate(y_full = change_p_i_full * y)

Now we compute the full endowment general equilibrium of aggregate expenditures (e_full and e_full_r):

ch2_application1 <- ch2_application1 %>% 
  mutate(e_full = change_p_j_full * e * (exporter == importer)) %>% 
  group_by(importer) %>% 
  mutate(e_full = max(e_full, na.rm = TRUE)) %>% 
  ungroup() %>% 
  mutate(
    e_full_r = e_full * (importer == "0-DEU"),
    e_full_r = max(e_full_r)
  )

With the aggregate expenditure we proceed to obtain the full endowment general equilibrium of the outward multilateral resistance (omr_full) and inward multilateral resistance (imr_full):

ch2_application1 <- ch2_application1 %>% 
  mutate(
    omr_full = y_full * e_r_cfl_0 / exp(fe_exporter_cfl_0),
    imr_full = e_full / (exp(fe_importer_cfl_0) * e_full_r)
  )

Finally we proceed to compute the full endowment general equilibrium of trade (xi_full):

ch2_application1 <- ch2_application1 %>% 
  mutate(x_full = (y_full * e_full * tij_cfl) / (imr_full * omr_full)) %>% 
  group_by(exporter) %>% 
  mutate(xi_full = sum(x_full * (importer != exporter))) %>% 
  ungroup()

This is our dataset at this point:

ch2_application1
## # A tibble: 4,761 x 64
##    exporter importer pair_id  year   trade   dist  cntg  lang  clny log_dist
##    <chr>    <chr>      <int> <int>   <dbl>  <dbl> <int> <int> <int>    <dbl>
##  1 ARG      ARG        12339  2006 32313.    534.     0     0     0     6.28
##  2 ARG      AUS            1  2006   108.  12045.     0     0     0     9.40
##  3 ARG      AUT            2  2006    25.2 11751.     0     0     0     9.37
##  4 ARG      BEL            4  2006   189.  11305.     0     0     0     9.33
##  5 ARG      BGR            3  2006    17.0 12116.     0     0     0     9.40
##  6 ARG      BOL            6  2006   392.   1866.     1     1     0     7.53
##  7 ARG      BRA            8  2006  6198.   2392.     1     0     0     7.78
##  8 ARG      CAN           10  2006   339.   9391.     0     0     0     9.15
##  9 ARG      CHE           12  2006   278.  11233.     0     0     0     9.33
## 10 ARG      CHL           14  2006  2456.   1157.     1     1     0     7.05
## # … with 4,751 more rows, and 54 more variables: intl <dbl>, y <dbl>, e <dbl>,
## #   e_r <dbl>, fe_exporter_bln <dbl>, fe_importer_bln <dbl>, tij_bln <dbl>,
## #   omr_bln <dbl>, imr_bln <dbl>, tradehat_bln <dbl>, xi_bln <dbl>,
## #   tij_cfl <dbl>, log_tij_cfl <dbl>, fe_exporter_cfl <dbl>,
## #   fe_importer_cfl <dbl>, omr_cfl <dbl>, imr_cfl <dbl>, tradehat_cfl <dbl>,
## #   xi_cfl <dbl>, change_tij <dbl>, phi <dbl>, change_p_i <dbl>,
## #   change_p_j <dbl>, trade_cfl <dbl>, omr_cfl_0 <dbl>, imr_cfl_0 <dbl>,
## #   change_imr_full_0 <dbl>, change_omr_full_0 <dbl>, change_p_i_0 <dbl>,
## #   change_p_j_0 <dbl>, fe_exporter_cfl_0 <dbl>, fe_importer_cfl_0 <dbl>,
## #   tradehat_0 <dbl>, e_r_cfl_0 <dbl>, trade_1 <dbl>, tradehat_1 <dbl>,
## #   y_cfl_1 <dbl>, e_cfl_1 <dbl>, e_r_cfl_1 <dbl>, change_p_i_1 <dbl>,
## #   change_p_j_1 <dbl>, omr_cfl_1 <dbl>, imr_cfl_1 <dbl>,
## #   change_omr_full_1 <dbl>, change_imr_full_1 <dbl>, change_p_i_full <dbl>,
## #   change_p_j_full <dbl>, y_full <dbl>, e_full <dbl>, e_full_r <dbl>,
## #   omr_full <dbl>, imr_full <dbl>, x_full <dbl>, xi_full <dbl>

3.1.5 Step IV: Collect, construct, and report indexes of interest

Box #1 from page 108 in Yotov et al. (2016) consists in to construct the percentage change of the general equilibrium indexes. The steps are direct, we need to compute the change in full endowment general equilibrium factory-gate price on export side (change_price_full), the change in conditional and full general equilibrium outward multilateral resistances (change_omr_*), and the change in conditional and full general equilibrium international trade (change_xi_*):

ch2_application1 <- ch2_application1 %>% 
  mutate(
    change_price_full = (change_p_i_full - 1) * 100,
    change_omr_cfl = (omr_cfl^(1 / (1 - sigma)) / omr_bln^(1 / (1 - sigma)) - 1) * 100,
    change_omr_full = (omr_full^(1 / (1 - sigma)) / omr_bln^(1 / (1 - sigma)) - 1) * 100,
    change_xi_cfl = (xi_cfl / xi_bln  - 1) * 100,
    change_xi_full = (xi_full / xi_bln - 1) * 100
  )

In addition to this, we need to something very similar for importers, in order to be able to recreate figure 7 later:

ch2_application1 <- ch2_application1 %>% 
  mutate(
    change_imr_full = -(imr_full^(1 / (1 - sigma)) / imr_bln^(1 / (1 - sigma)) - 1) * 100,
    rgdp = ((y_full / imr_full^(1 / (1 - sigma))) / (y / imr_bln^(1 / (1 - sigma))) - 1) * 100
  )

Now we are almost ready, except for some filters and re-scaling variables to plot:

ch2_application1
## # A tibble: 4,761 x 71
##    exporter importer pair_id  year   trade   dist  cntg  lang  clny log_dist
##    <chr>    <chr>      <int> <int>   <dbl>  <dbl> <int> <int> <int>    <dbl>
##  1 ARG      ARG        12339  2006 32313.    534.     0     0     0     6.28
##  2 ARG      AUS            1  2006   108.  12045.     0     0     0     9.40
##  3 ARG      AUT            2  2006    25.2 11751.     0     0     0     9.37
##  4 ARG      BEL            4  2006   189.  11305.     0     0     0     9.33
##  5 ARG      BGR            3  2006    17.0 12116.     0     0     0     9.40
##  6 ARG      BOL            6  2006   392.   1866.     1     1     0     7.53
##  7 ARG      BRA            8  2006  6198.   2392.     1     0     0     7.78
##  8 ARG      CAN           10  2006   339.   9391.     0     0     0     9.15
##  9 ARG      CHE           12  2006   278.  11233.     0     0     0     9.33
## 10 ARG      CHL           14  2006  2456.   1157.     1     1     0     7.05
## # … with 4,751 more rows, and 61 more variables: intl <dbl>, y <dbl>, e <dbl>,
## #   e_r <dbl>, fe_exporter_bln <dbl>, fe_importer_bln <dbl>, tij_bln <dbl>,
## #   omr_bln <dbl>, imr_bln <dbl>, tradehat_bln <dbl>, xi_bln <dbl>,
## #   tij_cfl <dbl>, log_tij_cfl <dbl>, fe_exporter_cfl <dbl>,
## #   fe_importer_cfl <dbl>, omr_cfl <dbl>, imr_cfl <dbl>, tradehat_cfl <dbl>,
## #   xi_cfl <dbl>, change_tij <dbl>, phi <dbl>, change_p_i <dbl>,
## #   change_p_j <dbl>, trade_cfl <dbl>, omr_cfl_0 <dbl>, imr_cfl_0 <dbl>,
## #   change_imr_full_0 <dbl>, change_omr_full_0 <dbl>, change_p_i_0 <dbl>,
## #   change_p_j_0 <dbl>, fe_exporter_cfl_0 <dbl>, fe_importer_cfl_0 <dbl>,
## #   tradehat_0 <dbl>, e_r_cfl_0 <dbl>, trade_1 <dbl>, tradehat_1 <dbl>,
## #   y_cfl_1 <dbl>, e_cfl_1 <dbl>, e_r_cfl_1 <dbl>, change_p_i_1 <dbl>,
## #   change_p_j_1 <dbl>, omr_cfl_1 <dbl>, imr_cfl_1 <dbl>,
## #   change_omr_full_1 <dbl>, change_imr_full_1 <dbl>, change_p_i_full <dbl>,
## #   change_p_j_full <dbl>, y_full <dbl>, e_full <dbl>, e_full_r <dbl>,
## #   omr_full <dbl>, imr_full <dbl>, x_full <dbl>, xi_full <dbl>,
## #   change_price_full <dbl>, change_omr_cfl <dbl>, change_omr_full <dbl>,
## #   change_xi_cfl <dbl>, change_xi_full <dbl>, change_imr_full <dbl>,
## #   rgdp <dbl>

3.1.6 Figures replication

With all of the steps above, we are ready to create the plots from page 110. in Yotov et al. (2016).

Figure 6 removes the observations where both the importer and the exporter are different, this can be seen in the original Stata code provided with the book.

We need to filter rows and to obtain log(y):

ch2_application1 <- ch2_application1 %>% 
  filter(exporter == importer) %>% 
  select(exporter, importer, y, change_xi_cfl, change_xi_full, rgdp, 
         change_price_full, change_imr_full) %>% 
  mutate(log_y = log(y))

In addition, the original code removes Hong Kong for visualization scale purposes:

ggplot(data = ch2_application1 %>% 
         filter(exporter != "HKG")) +
  geom_point(aes(x = log_y, y = change_xi_cfl, color = "1")) + 
  geom_point(aes(x = log_y, y = change_xi_full, color = "2")) +
  labs(
    x = "Log value of output",
    y = "Percent change of exports",
    title = "Figure 6: Effects of abolishing international borders on exports",
    caption = "Source: Authors' calculations",
    color = ""
  ) +
  theme_minimal() +
  theme(legend.position = "bottom") +
  scale_color_manual(
    labels = c(
      "Conditional general equilibrium",
      "Full endowment general equilibrium"
    ),
    values = c("#b6b8dd","#232958")
  )

To create figure 7, we proceed in the same way as we did with figure 6:

ggplot(data = ch2_application1 %>% filter(exporter != "HKG")) +
  geom_point(aes(x = log_y, y = change_imr_full, color = "1")) + 
  geom_point(aes(x = log_y, y = change_price_full, color = "2")) +
  geom_point(aes(x = log_y, y = rgdp, color = "3")) +
  labs(
    x = "Log value of output",
    y = "Percent changes",
    title = "Figure 7: Effects of abolishing international borders on real GDP",
    caption = "Note: The inward multilateral resistances have been reformulated by multiplying their value by minus one.\nSource: Authors' calculations",
    color = ""
  ) +
  theme_minimal() +
  theme(legend.position = "bottom") +
  scale_color_manual(
    labels = c(
      "-(inward multilateral resistances)",
      "Factory-gate price",
      "Real GDP"
    ),
    values = c("#3bade3", "#b6b8dd", "#232958")
  )

3.2 Impact of regional trade agreements

3.2.1 Initial data

As in the previous application, we shall proceed by alternating both data transforming and regressions. Before doing what is shown in box #1 from page 112 in Yotov et al. (2016), we need to convert “DEU” in both exporter and importer columns to “0-DEU,” just as in the last section, but we shall keep panel dimension of the dataset in order to identify the effects of RTAs and to comprehensively capture the impact of all time-invariant trade costs with the use of pair fixed effects:

ch2_application2 <- yotov_data("ch2_application2") %>% 
  filter(year %in% seq(1986, 2006, 4)) %>% 
  mutate(
    log_dist = log(dist),
    intl = ifelse(exporter != importer, 1, 0),
    exporter = ifelse(exporter == "DEU", "0-DEU", exporter),
    importer = ifelse(importer == "DEU", "0-DEU", importer)
  ) %>% 
  
  # Create Yit
  group_by(exporter, year) %>% 
  mutate(y = sum(trade)) %>% 
  
  # Create Eit
  group_by(importer, year) %>% 
  mutate(e = sum(trade)) %>% 
  
  # Create Er
  group_by(year) %>% 
  mutate(e_r = max(ifelse(importer == "0-DEU", e, NA), na.rm = T))

Because of the panel dimension, we proceed as we did in the previous chapter, by creating columns to combine exporter/importer and year (exp_year and imp_year) for the fixed effects, and a pairing variable pair_id_2:

ch2_application2 <- ch2_application2 %>% 
  mutate(
    exp_year = paste0(exporter, year),
    imp_year = paste0(importer, year),
    pair_id_2 = ifelse(exporter == importer, "0-intra", pair_id)
  )

In addition, we need to create the variable sum_trade to filter the cases where the sum by pair_id is zero:

ch2_application2 <- ch2_application2 %>% 
  group_by(pair_id) %>% 
  mutate(sum_trade = sum(trade)) %>% 
  ungroup()

3.2.2 Step 1: Solve the baseline gravity model

3.2.2.1 Stage 1: Obtain the estimates of pair fixed effects and the effects of RTAs

With the steps done before, it is straightforward to obtain the PPML regression shown in box #1 from page 112 in Yotov et al. (2016):

fit_baseline_app2 <- glm(
  trade ~ 0 + rta + exp_year + imp_year + pair_id_2,
  family = quasipoisson(link = "log"),
  data = filter(ch2_application2, sum_trade > 0)
)

With the estimated model, we can construct the variables for export (exp_year), import (imp_year) and pair (pair_id_2) fixed effects:

ch2_application2 <- ch2_application2 %>%
  left_join(
    yotov_fixed_effects(fit_baseline_app2),
    by = c("exp_year", "imp_year", "pair_id_2")
  )

3.2.2.2 Stage 2: Regress the estimates of pair fixed effects on gravity variables and country fixed effects

Box #1 from page 113 in Yotov et al. (2016) can be divided in smaller chunks.

We start by filtering to keep the observation from 1994, and then we compute the trade costs (tij_bar and tijn_bln) from the fe_pair_id_2 fixed effects and the estimated RTA coefficient:

ch2_application2 <- ch2_application2 %>% 
  mutate(
    tij_bar = exp(fe_pair_id_2),
    tij_bln = exp(fe_pair_id_2 + fit_baseline_app2$coefficients["rta"] * rta)
  )

Now we need to create a table for the year 1994 which will be used to predict the trade costs for the observations with zero trade flows. The reason to create a sub-table, instead of filtering observations in the regression function, is that is eases posterior work to predict the costs.

To start the cost estimation we start with:

ch2_application2_1994 <- ch2_application2 %>% 
  filter(year == 1994, exporter != importer)

Now, unlike the book, which duplicates tij_bar by creating tij, we can fit a regression to estimate the costs:

fit_costs_app2 <- glm(
  tij_bar ~ 0 + log_dist + cntg + lang + clny + exporter + importer,
  family = quasipoisson(link = "log"),
  data = ch2_application2_1994
)

With the regression we add the fitted values to the sub-table:

ch2_application2_1994 <- ch2_application2_1994 %>% 
  mutate(tij_no_rta = predict(fit_costs_app2, ch2_application2_1994, "response")) %>% 
  select(exporter, importer, tij_no_rta)

The final step is to keep the observations for the year 1994 in the original table and replace the missing costs with the predicted values:

ch2_application2 <- ch2_application2 %>% 
  filter(year == 1994) %>% 
  left_join(ch2_application2_1994, by = c("exporter", "importer")) %>% 
  mutate(
    tij_bar = ifelse(is.na(tij_bar), tij_no_rta, tij_bar),
    tij_bln = ifelse(is.na(tij_bln), tij_bar * exp(fit_baseline_app2$coefficients["rta"] * rta), tij_bln)
  ) %>% 
  select(-tij_no_rta) %>% 
  mutate(log_tij_bln = log(tij_bln))

Box #2 from page 113 in Yotov et al. (2016) is more straightforward.

The first part to complete Box #2 consists in solving the constrained baseline gravity model:

fit_constrained_app2 <- glm(
  trade ~ 0 + exporter + importer + offset(log_tij_bln),
  family = quasipoisson(link = "log"),
  data = ch2_application2
)

With the fitted model we can add the prediction and the xi_bln column:

ch2_application2 <- ch2_application2 %>% 
  mutate(tradehat_bln = predict(fit_constrained_app2, ch2_application2, "response")) %>% 
  group_by(exporter) %>% 
  mutate(xi_bln = sum(tradehat_bln * (exporter != importer))) %>% 
  ungroup()

The book specifies that all other baseline indexes of interest can be obtained by applying the exact same procedure as described in the previous application. Here we’ll obtain the multilateral resistances terms (omr_bln and imr_bln) by adding the fixed effects from the constrained model to the data:

ch2_application2 <- ch2_application2 %>% 
  left_join(yotov_fixed_effects(fit_constrained_app2), by = c("exporter","importer"))

ch2_application2 <- ch2_application2 %>% 
  mutate(
    omr_bln = y * e_r/ exp(fe_exporter),
    imr_bln = e / (exp(fe_importer) * e_r)
  )

3.2.3 Step II: Define a counterfactual scenario

Box #1 from page 114 in Yotov et al. (2016) is direct and consists in replacing the RTA values by zero if the pairs of countries are NAFTA members:

nafta <- c("MEX", "USA", "CAN")

ch2_application2 <- ch2_application2 %>% 
  mutate(
    rta_no_nafta = ifelse(exporter %in% nafta & importer %in% nafta, 0, rta),
    tij_cfl = tij_bar * exp(fit_baseline_app2$coefficients["rta"] * rta_no_nafta),
    log_tij_cfl = log(tij_cfl)
  )

3.2.4 Step III: Solve the counterfactual model

The exact same procedure from the previous section applies to obtain the conditional general equilibrium effects and then to compute the full endowment general equilibrium effects.

We start by fitting a counterfactual model:

fit_counterfactual_app2 <- glm(
  trade ~ 0 + exporter + importer + offset(log_tij_cfl),
  family = quasipoisson(link = "log"),
  data = ch2_application2
)

With the fitted model we add the fixed effects:

ch2_application2 <- ch2_application2 %>% 
  left_join(yotov_fixed_effects(fit_counterfactual_app2), by = c("exporter","importer")) %>% 
  rename(
    fe_exporter_bln = fe_exporter.x,
    fe_exporter_cfl = fe_exporter.y,
    fe_importer_bln = fe_importer.x,
    fe_importer_cfl = fe_importer.y
  )

As we did in stage 2, we compute the multilateral resistance terms:

ch2_application2 <- ch2_application2 %>% 
  mutate(
    omr_cfl = y * e_r / exp(fe_exporter_cfl),
    imr_cfl = e / (exp(fe_importer_cfl) * e_r)
  )

Up to this point we are ready to compute the conditional general equilibrium effects of trade:

ch2_application2 <- ch2_application2 %>% 
  mutate(tradehat_cfl = predict(fit_counterfactual_app2, ch2_application2, "response")) %>% 
  group_by(exporter) %>% 
  mutate(xi_cfl = sum(tradehat_cfl * (exporter != importer))) %>% 
  ungroup()

We are going to compute the full endowment general equilibrium effects. We start by repeating the steps to obtain change_tij and phi from the last section:

# set the criteria of convergence
# taken from the literature (see the Stata code)
sigma <- 7

ch2_application2 <- ch2_application2 %>% 
  mutate(
    change_tij = tij_cfl / tij_bln,
    phi = ifelse(importer == exporter, e / y, 0)
  ) %>% 
  group_by(exporter) %>% 
  mutate(phi = max(phi)) %>% 
  ungroup()

Now we compute change_p_i, change_p_j and trade_cfl. Again, this is just a repetition of the previous steps with some adaptation:

ch2_application2 <- ch2_application2 %>% 
  group_by(exporter) %>% 
  mutate(change_p_i = ((exp(fe_exporter_cfl) / e_r) / (exp(fe_exporter_bln) / e_r))^(1 /(1 - sigma))) %>% 
  ungroup() %>% 
  
  group_by(importer) %>% 
  mutate(
    change_p_j = ifelse(importer == exporter, change_p_i, 0),
    change_p_j = max(change_p_j)
  ) %>% 
  ungroup()

ch2_application2 <- ch2_application2 %>% 
  mutate(trade_cfl = tradehat_cfl * change_p_i * change_p_j)

Then we need a while() loop, but before we need to duplicate some columns under new names for the loop operations:

ch2_application2 <- ch2_application2 %>% 
  mutate(
    omr_cfl_0 = omr_cfl,
    imr_cfl_0 = imr_cfl,
    change_imr_full_0 = 1,
    change_omr_full_0 = 1,
    change_p_i_0 = change_p_i,
    change_p_j_0 = change_p_j,
    fe_exporter_cfl_0 = fe_exporter_cfl,
    fe_importer_cfl_0 = fe_importer_cfl,
    tradehat_0 = tradehat_cfl,
    e_r_cfl_0 = e_r
  )

Now we run the loop, where the \(N^{th}\) step depends on the \(1,2,\ldots,N-1\) previous steps as in the previous section:

# set parameters
max_dif <- 1
sd_dif <- 1
change_price_i_old <- 0

i2 <- 1
while(sd_dif > 1e-3 | max_dif > 1e-3) {
  ch2_application2 <- ch2_application2 %>% 
    mutate(trade_1 = tradehat_0 * change_p_i_0 * change_p_j_0 / (change_omr_full_0 * change_imr_full_0))
  
  # repeat the counterfactual model
  fit_counterfactual_app2_2 <- glm(
    trade_1 ~ 0 + exporter + importer + offset(log_tij_cfl),
    family = quasipoisson(link = "log"),
    data = ch2_application2
  )
  
  ch2_application2 <- ch2_application2 %>% 
    left_join(
      yotov_fixed_effects(fit_counterfactual_app2_2),
      by = c("exporter", "importer")
    )
  
  # compute the conditional general equilibrium effects of trade
  ch2_application2 <- ch2_application2 %>% 
    mutate(tradehat_1 = predict(fit_counterfactual_app2_2, ch2_application2, "response")) %>% 
    group_by(exporter) %>% 
    mutate(y_cfl_1 = sum(tradehat_1)) %>% 
    ungroup() %>% 
    
    mutate(e_cfl_1 = ifelse(importer == exporter, phi * y_cfl_1, 0)) %>% 
    group_by(importer) %>% 
    mutate(e_cfl_1 = max(e_cfl_1)) %>% 
    ungroup() %>% 
    
    mutate(
      e_r_cfl_1 = ifelse(importer == "0-DEU", e_cfl_1, 0),
      e_r_cfl_1 = max(e_r_cfl_1)
    )
  
  # compute the change in prices for exporters and importers
  ch2_application2 <- ch2_application2 %>% 
    mutate(change_p_i_1 = ((exp(fe_exporter) / e_r_cfl_1) / 
      (exp(fe_exporter_cfl_0) / e_r_cfl_0))^(1 / (1 - sigma)))
  
  # compute the change in prices for exporters and importers
  ch2_application2 <- ch2_application2 %>% 
    group_by(importer) %>% 
    mutate(
      change_p_j_1 = ifelse(importer == exporter, change_p_i_1, 0),
      change_p_j_1 = max(change_p_j_1)
    ) %>% 
    ungroup()
  
  # compute both outward and inward multilateral resistance
  ch2_application2 <- ch2_application2 %>% 
    mutate(
      omr_cfl_1 = (y_cfl_1 * e_r_cfl_1) / exp(fe_exporter),
      imr_cfl_1 = e_cfl_1 / (exp(fe_importer) * e_r_cfl_1)
    )
  
  # update the differences
  max_dif <- abs(max(ch2_application2$change_p_i_0 - change_price_i_old))
  sd_dif <- sd(ch2_application2$change_p_i_0 - change_price_i_old)
  change_price_i_old <- ch2_application2$change_p_i_0
  
  # compute changes in outward and inward multilateral resistance
  ch2_application2 <- ch2_application2 %>% 
    mutate(
      change_omr_full_1 = omr_cfl_1 / omr_cfl_0,
      change_imr_full_1 = imr_cfl_1 / imr_cfl_0,
      omr_cfl_0 = omr_cfl_1,
      imr_cfl_0 = imr_cfl_1,
      change_omr_full_0 = change_omr_full_1,
      change_imr_full_0 = change_imr_full_1,
      change_p_i_0 = change_p_i_1,
      change_p_j_0 = change_p_j_1,
      fe_exporter_cfl_0 = fe_exporter,
      fe_importer_cfl_0 = fe_importer,
      tradehat_0 = tradehat_1,
      e_r_cfl_0 = e_r_cfl_1
    ) %>% 
    select(-fe_exporter, -fe_importer)
  
  i2 <- i2 + 1
}

The last loop allows us to obtain change_p_i_full, change_p_j_full and y_full:

ch2_application2 <- ch2_application2 %>%
  mutate(
    change_p_i_full = ((exp(fe_exporter_cfl_0) / e_r_cfl_0) /
                         (exp(fe_exporter_bln) / e_r))^(1 / (1 - sigma)),
    change_p_j_full = change_p_i_full * (exporter == importer)
  ) %>%
  group_by(importer) %>%
  mutate(change_p_j_full = max(change_p_j_full)) %>%
  ungroup() %>%
  mutate(y_full = change_p_i_full * y)

Now we compute e_full and e_full_r:

ch2_application2 <- ch2_application2 %>%
  mutate(e_full = change_p_j_full * e * (exporter == importer)) %>%
  group_by(importer) %>%
  mutate(e_full = max(e_full, na.rm = TRUE)) %>%
  ungroup() %>%
  mutate(
    e_full_r = e_full * (importer == "0-DEU"),
    e_full_r = max(e_full_r)
  )

We also need omr_full and imr_full. This part of the code needs attention because the IMR full term is computed in a different way compared to the previous section, please see the script RTAsEffects.do. The way to replicate the original Stata code is:

ch2_application2 <- ch2_application2 %>%
  mutate(
    omr_full = y_full * e_r_cfl_0 / exp(fe_exporter_cfl_0),
    imr_full = e_cfl_1 / (exp(fe_importer_cfl_0) * e_r_cfl_0)
  )

To complete this step, we compute xi_full:

ch2_application2 <- ch2_application2 %>%
  mutate(x_full = (y_full * e_full * tij_cfl) / (imr_full * omr_full)) %>%
  group_by(exporter) %>%
  mutate(xi_full = sum(x_full * (importer != exporter))) %>%
  ungroup()

3.2.5 Step IV: Collect, construct, and report indexes of interest

The goal of this step is to reproduce the table from page 116 in Yotov et al. (2016).

To ease the task of creating the table from the book, we divide between exporter and importer indexes:

exporter_indexes <- ch2_application2 %>% 
  select(
    exporter, starts_with("omr_"), change_p_i_full,
    starts_with("xi_"), y, y_full
  ) %>% 
  distinct() %>% 
  mutate(exporter = ifelse(exporter == "0-DEU", "DEU", exporter)) %>% 
  arrange(exporter) %>% 
  mutate(
    change_p_i_full = (1 - change_p_i_full) * 100,
    change_omr_cfl = ((omr_bln / omr_cfl)^(1 / (1-sigma)) - 1) * 100,
    change_omr_full = ((omr_bln / omr_full)^(1 / (1-sigma)) - 1) * 100,
    change_xi_cfl = (xi_bln / xi_cfl - 1) * 100,
    change_xi_full = (xi_bln / xi_full - 1) * 100
  ) %>% 
 select(exporter, starts_with("change"), starts_with("y"))
  
importer_indexes <- ch2_application2 %>% 
  select(importer, imr_bln, imr_cfl, imr_full) %>% 
  ungroup()  %>% 
  distinct() %>% 
  mutate(importer = ifelse(importer == "0-DEU", "DEU", importer)) %>% 
  arrange(importer) %>% 
  mutate(
    change_imr_cfl = ((imr_bln / imr_cfl)^(1 / (1 - sigma)) - 1) * 100,
    change_imr_full = ((imr_bln / imr_full)^(1 / (1 - sigma)) - 1) * 100
  )

Finally, we can replicate the table that we wanted to:

indexes_final <- exporter_indexes %>% 
  left_join(importer_indexes, by = c("exporter" = "importer")) %>% 
  mutate(
    rgdp_bln = y / (imr_bln^(1 / (1 - sigma))),
    rgdp_full = y_full / (imr_full^(1 / (1 - sigma))),
    change_rgdp_full = (rgdp_bln / rgdp_full - 1) * 100
  ) %>% 
  select(exporter, change_xi_cfl, change_xi_full, 
         change_rgdp_full, change_imr_full, change_omr_full, change_p_i_full)

indexes_final %>% 
  mutate_if(is.numeric, function(x) round(x, 2))
## # A tibble: 69 x 7
##    exporter change_xi_cfl change_xi_full change_rgdp_full change_imr_full
##    <chr>            <dbl>          <dbl>            <dbl>           <dbl>
##  1 ARG              -0.66          -0.69            -0.01            0.01
##  2 AUS              -0.49          -0.51            -0.01            0.02
##  3 AUT              -0.07          -0.1             -0.01            0   
##  4 BEL              -0.09          -0.12             0               0   
##  5 BGR              -0.05          -0.07             0               0   
##  6 BOL              -0.47          -0.48            -0.02            0.03
##  7 BRA              -0.65          -0.69            -0.01           -0.01
##  8 CAN              35.0           37.5              3.4            -1.48
##  9 CHE              -0.14          -0.17            -0.01            0   
## 10 CHL              -0.81          -0.83            -0.03            0.01
## # … with 59 more rows, and 2 more variables: change_omr_full <dbl>,
## #   change_p_i_full <dbl>

References

Yotov, Yoto V, Roberta Piermartini, José-Antonio Monteiro, and Mario Larch. 2016. An Advanced Guide to Trade Policy Analysis: The Structural Gravity Model. World Trade Organization Geneva.