01-traditional-gravity.qmd

---
bibliography: 00-references.bib
---

# Partial equilibrium trade policy analysis with structural gravity

## Traditional Gravity Estimates

### Preparing the data

If the reader has never used R before, please check chapters 1 to 25 from @wickham2016r.

If the reader has only fitted a few regressions in R, without much practice on transforming and cleaning data before, please check chapters 5 and 18 from @wickham2016r.

Please see the note from page 42 in @yotov2016advanced. It is a really important note, which tells us that we need to:

1. Filter observations for a range of years (1986, 1990, 1994, 1998, 2002 and 2006)
2. Transform some variables to logarithm scale (trade and dist) and create new variables from those in the original dataset
3. Remove cases where both the exporter and the importer are the same
4. Drop observations where the trade flow is zero

Unlike @yotov2016advanced, here we shall use a single dataset for all the applications and subset its columns depending on what we need. This decision kept the *tradepolicy* R package as light as possible.

Before conducting any data filtering or regression, we need to load the required packages.

```{r ch1_app_1_packages}
#| message: false
#| warning: false

# dataset and summary functions
library(tradepolicy)

# data transformation
library(dplyr)
library(tidyr)

# regression
library(fixest)
```

Step 1, including subsetting columns for this application, is straightforward.

```{r ch1_app_1_data_1}
ch1_application1 <- agtpa_applications %>%
  select(exporter, importer, pair_id, year, trade, dist, cntg, lang, clny) %>%
  filter(year %in% seq(1986, 2006, 4))
```

For step 2, this can be divided in parts, starting with the log transformation of trade and distance.

```{r ch1_app_1_data_2}
ch1_application1 <- ch1_application1 %>%
  mutate(
    log_trade = log(trade),
    log_dist = log(dist)
  )
```

Continuing step 2, we can now create the variables $Y_{i,t}$ and $E_{i,t}$ that appear on the OLS model equation in the book.

```{r ch1_app_1_data_3}
ch1_application1 <- ch1_application1 %>%
  # Create Yit
  group_by(exporter, year) %>%
  mutate(
    y = sum(trade),
    log_y = log(y)
  ) %>%
  # Create Eit
  group_by(importer, year) %>%
  mutate(
    e = sum(trade),
    log_e = log(e)
  )
```

The OLS model with remoteness index needs both exporter and importer index, which grouping variables can create. We divide it into sub-steps: Replicate the computation of total exports, then the remoteness index for exporters, and finally the total imports with the corresponding remoteness index for importers.

```{r ch1_app_1_data_4}
ch1_application1 <- ch1_application1 %>%
  # Replicate total_e
  group_by(exporter, year) %>%
  mutate(total_e = sum(e)) %>%
  group_by(year) %>%
  mutate(total_e = max(total_e)) %>%
  # Replicate rem_exp
  group_by(exporter, year) %>%
  mutate(
    remoteness_exp = sum(dist * total_e / e),
    log_remoteness_exp = log(remoteness_exp)
  ) %>%
  # Replicate total_y
  group_by(importer, year) %>%
  mutate(total_y = sum(y)) %>%
  group_by(year) %>%
  mutate(total_y = max(total_y)) %>%
  # Replicate rem_imp
  group_by(importer, year) %>%
  mutate(
    remoteness_imp = sum(dist / (y / total_y)),
    log_remoteness_imp = log(remoteness_imp)
  )
```

To create the variables for the OLS with Fixed Effects Model, we followed box #1 on page 44 from @yotov2016advanced. We combine both exporter and importer variables with the year to create the fixed effects variables.

```{r ch1_app_1_data_5}
ch1_application1 <- ch1_application1 %>%
  # This merges the columns exporter/importer with year
  mutate(
    exp_year = paste0(exporter, year),
    imp_year = paste0(importer, year)
  )
```

The addition of exporter/importer time fixed effects concludes step 2, and now we need to perform step 3.

```{r ch1_app_1_data_6}
ch1_application1 <- ch1_application1 %>%
  filter(exporter != importer)
```

Some cases require conducting step 4, and we will be explicit about it when needed.

### OLS estimation ignoring multilateral resistance terms

The general equation for this model is
$$
\begin{align}
\log X_{ij,t} =& \:\beta_0 + \beta_1 DIST_{i,j} + \beta_2 CNTG_{i,j} + \beta_3 LANG_{i,j} + \beta_4 CLNY_{i,j} + \beta_5 \log Y_{i,t} +\\
\text{ }& \:\beta_6 \log E_{j,t} + \varepsilon_{ij,t}.
\end{align}
$$

Please see page 41 in @yotov2016advanced for full detail of each variable.

The model for this case is straightforward, and in this case, we need to apply step 4 from the previous section to drop cases where the trade is zero.

```{r ch1_app_1_ols_1}
fit_ols <- feols(
  log_trade ~ log_dist + cntg + lang + clny + log_y + log_e,
  data = filter(ch1_application1, trade > 0)
)

summary(fit_ols)
```

The employed function, `feols()`, does not carry a copy of its training data by default besides providing faster fitting for models with fixed effects. This is not the case in base R, where `glm()` outputs include this data, increasing the model's size, but this does not affect the model's predictions and can be changed as the user needs it [@trimmingfat].

The model is almost ready. We only need to stick to the methodology from @yotov2016advanced and cluster the standard errors by country pair (see the note on page 42, it is imperative).

```{r ch1_app_1_ols_1_2}
fit_ols <- feols(
  log_trade ~ log_dist + cntg + lang + clny + log_y + log_e,
  data = filter(ch1_application1, trade > 0),
  cluster = ~pair_id
)

summary(fit_ols)
```

The *tradepolicy* package provides functions to provide more informative summaries. Please read the documentation of the package and look for the `tp_summary_app_1()` function, it summarises the model in the exact way as reported in the book by providing:

* Clustered standard errors.
* Number of observations.
* $R^2$ (if applicable).
* Presence (or absence) of exporter and exporter-time fixed effects.
* RESET test p-value.

These statistical results are returned as a list to keep it simple, which we can see for the model in the same format as reported in the book.

```{r ch1_app_1_ols_2}
tp_summary_app_1(
  formula = log_trade ~ log_dist + cntg + lang + clny + log_y + log_e,
  data = filter(ch1_application1, trade > 0),
  method = "ols"
)
```

Please notice that the summary hides the exporter/importer fixed effects.

### OLS estimation controlling for multilateral resistance terms with remote indexes

The remoteness model adds variables to the OLS model. The general equation for this model is
$$
\begin{align}
\log X_{ij,t} =& \:\beta_0 + \beta_1 DIST_{i,j} + \beta_2 CNTG_{i,j} + \beta_3 LANG_{i,j} + \beta_4 CLNY_{i,j} + \beta_5 \log Y_{i,t} +\\
\text{ }& \beta_6 \log E_{j,t} + \beta_7 \log(REM\_EXP_i,t) + \beta_8 \log(REM\_IMP_i,t) + \varepsilon_{ij,t}.
\end{align}
$$

In the equation above $REM\_EXP$ and $REM\_IMP$ are defined as
$$
\begin{align}
\log(REM\_EXP_{i,t}) &= \log \left( \sum_j \frac{DIST_{i,j}}{E_{j,t} / Y_t} \right) \text{ and }\\
\log(REM\_IMP_{j,t}) &= \log \left( \sum_i \frac{DIST_{i,j}}{Y_{i,t} / Y_t} \right).
\end{align}
$$

Please see page 43 in @yotov2016advanced for full detail of each variable.

Our approach follows box #1 on page 43 from @yotov2016advanced. Fitting the regression is straightforward. It is just about adding more regressors to what we did in the last section, and we can create a list with a summary for the model.

```{r ch1_app_1_ols_remoteness_1}
tp_summary_app_1(
  formula = log_trade ~ log_dist + cntg + lang + clny + log_y + log_e +
    log_remoteness_exp + log_remoteness_imp,
  data = filter(ch1_application1, trade > 0),
  method = "ols"
)
```

### OLS estimation controlling for multilateral resistance terms with fixed effects

The general equation for this model is
$$
\begin{align}
\log X_{ij,t} =& \: \beta_1 \log(DIST)_{i,j} + \beta_2 CNTG_{i,j} + \beta_3 LANG_{i,j} +\\
\text{ }& \:\beta_4 CLNY_{i,j} + \pi_{i,t} + \chi_{i,t} + \varepsilon_{ij,t}.
\end{align}
$$

Where the added terms, concerning the OLS model, are $\pi_{i,t}$ and $\chi_{i,t}$ that account for exporter-time and importer-time fixed effects, respectively. See page 44 in @yotov2016advanced for full detail of each variable.

We can quickly generate a list as we did with the previous models. The only difference to the previous models is that in this case that the variables to the right of the "|" operator are the fixed effects, which are treated differently by the *fixest* package, which is used internally by the *tradepolicy* package, for faster model fitting.

Please notice that the summaries intentionally do not show fixed effects, because there are cases where we have thousands of fixed effects.

```{r ch1_app_1_fe_1}
tp_summary_app_1(
  formula = log_trade ~ log_dist + cntg + lang + clny | exp_year + imp_year,
  data = filter(ch1_application1, trade > 0),
  method = "ols"
)
```

There is another difference when we compare `feols()` or `fepois()` against `glm()` in the presence of fixed effects, which we can explain with an example.

In the data used for the previous summary, we have $T$ years (1986, 1990, 1994, 1998, 2002 and 2006). We could be interested in filtering for a single exporter and a single importer to fit the fixed effects model
$$
\log X_{t} = \beta_1 \log(DIST)_{i,j} + \beta_2 CNTG_{i,j} + \beta_3 LANG_{i,j} + \beta_4 CLNY_{i,j} + \sum_{u} \beta_u FE_u + \varepsilon_{ij,t},
$$
where
$$
\begin{align*}
&u \in \{1986, 1990, 1994, 1998, 2002, 2006\} \text{ and}\cr
&FE_u = \begin{cases}
1 & \text {if } t = u \cr
0 & \text{otherwise}.
\end{cases}
\end{align*}
$$

When we use `feols()`, or any of the functions in the *fixest* package, a formula of the form $z \sim x_1 + x_2 \mid y$ will estimate the model described in the previous equations.

If we do the same in base R, with `glm()`, the equivalent formula would be of the form $z \sim 0 + x_1 + x_2 + y$, otherwise base R estimates a model of the form $y_t = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \sum_u y_{u}$, which means that not including a zero in the formula estimates an additional coefficient in the estimation which corresponds to a global intercept or the "grand mean".

On the $\beta_0$ constant, it is tough to interpret with many-way fixed effects, and it is advised against reporting it.

### PPML estimation controlling for multilateral resistance terms with fixed effects

The general equation for this model is
$$
\begin{align}
X_{ij,t} =& \:\exp\left[\beta_1 \log(DIST)_{i,j} + \beta_2 CNTG_{i,j} +\right.\\
\text{ }& \:\left.\beta_3 LANG_{i,j} + \beta_4 CLNY_{i,j} + \pi_{i,t} + \chi_{i,t}\right] \times \varepsilon_{ij,t}.
\end{align}
$$

The reason to compute this model, despite the lower speed compared to OLS, is that PPML is the only estimator perfectly consistent with the theoretical gravity model. By estimating with PPML, the fixed effects correspond precisely to the corresponding theoretical terms.

The data for this model is the same as for the fixed effects model, and one option in R is to use the `fepois()` function.

```{r ch1_app_1_ppml_1, eval = FALSE}
fit_ppml <- fepois(trade ~ log_dist + cntg + lang + clny | exp_year + imp_year,
  data = ch1_application1,
  cluster = ~pair_id
)
```

Exactly as it was mentioned for `feols()`, `fepois()` shares the same differences regarding `glm()` objects.

If the reader decides to run this model and print the summary, they will notice that it does not report an $R^2$ and displays an extensive fixed effect list. For PPML models, the $R^2$ needs to be computed afterwards as Kendall's correlation (or rank correlation) between the observed and predicted values. Please see @silva2006log for the details and the RESET test for PPML (GLM) models.

Beware that software such as Stata requires additional libraries such as *ppmlhdfe* to report a correct $R^2$ for the PPML model. What Stata shows is a reported pseudo-$R^2$. To construct a proper $R^2$ in R, `tp_summary_app_1()` takes the rank correlation between actual and predicted trade flows.

We can obtain a detailed list as in the previous examples.

```{r ch1_app_1_ppml_2}
tp_summary_app_1(
  formula = trade ~ log_dist + cntg + lang + clny | exp_year + imp_year,
  data = ch1_application1,
  method = "ppml"
)
```

## The "distance puzzle" resolved

### Preparing the data

Please see the note from page 47 in @yotov2016advanced. We need to proceed with similar steps as in the previous section.

The distance puzzle proposes the gravity specification
$$
\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 LANG_{i,j}\right]\times\\
\text{ }& \:\exp\left[\beta_4 CLNY_{i,j} + \beta_5 \log(DIST\_INTRA_{i,i})\right] \times \varepsilon_{ij,t}.
\end{align}
$$

The difference concerning the last section is that now we need to separate the distance variable into multiple columns that account for discrete-time effects. The $\beta_T$ terms of the equation reflect this. Perhaps the easiest option is to transform the year into a text column and then use the `pivot_wider()` function.

We need to remove cases where the exporter is the same as the importer and cases where trade is zero for the OLS model. For the PPML models, we need to mark rows where the exporter and the importer are the same, and we need to create the same country column, which is also required to transform the distance variables as shown in box #1 in page 48 from @yotov2016advanced.

In order to avoid creating two very similar datasets, we shall create one dataset to cover both OLS and PPML.

```{r ch1_app_2_data_1}
ch1_application2 <- agtpa_applications %>%
  select(exporter, importer, pair_id, year, trade, dist, cntg, lang, clny) %>%
  # this filter covers both OLS and PPML
  filter(year %in% seq(1986, 2006, 4)) %>%
  mutate(
    # variables for both OLS and PPML
    exp_year = paste0(exporter, year),
    imp_year = paste0(importer, year),
    year = paste0("log_dist_", year),
    log_trade = log(trade),
    log_dist = log(dist),
    smctry = ifelse(importer != exporter, 0, 1),

    # PPML specific variables
    log_dist_intra = log_dist * smctry,
    intra_pair = ifelse(exporter == importer, exporter, "inter")
  ) %>%
  pivot_wider(names_from = year, values_from = log_dist, values_fill = 0) %>%
  mutate(across(log_dist_1986:log_dist_2006, function(x) x * (1 - smctry)))
```

The `across()` function is a shortcut to avoid repetition, as in the following example, we show it for reference without computation.

```{r ch1_app_2_data_3, eval = FALSE}
ch1_application2 %>%
  mutate(
    log_dist_1986 =  log_dist_1986 * (1 - smctry),
    log_dist_1990 =  log_dist_1990 * (1 - smctry),

    # repeat log_dist_T many_times for T = 1994, 1998, ...

    log_dist_2006 =  log_dist_2006 * (1 - smctry)
  )
```

Note that the OLS model shall require filtering when we specify the model because we skipped filtering the cases where trade is equal to zero and both the importer and the exporter are the same. Because the solution for the "distance puzzle" implies different transformations and filters for the OLS and PPML cases, one possibility is to filter in the same summary functions.

### OLS solution for the "distance puzzle"

The gravity specification, which includes $\pi_{i,t} + \chi_{i,t}$, means that we need to do something very similar to what we did in the last section.

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

```{r ch1_app_2_ols_1}
tp_summary_app_2(
  formula = log_trade ~ log_dist_1986 + log_dist_1990 + log_dist_1994 +
    log_dist_1998 + log_dist_2002 + log_dist_2006 + cntg + lang + clny |
    exp_year + imp_year,
  data = filter(ch1_application2, importer != exporter, trade > 0),
  method = "ols"
)
```

### PPML solution for the "distance puzzle"

This model is very similar to the one specified in the PPML section from the last section. We can directly fit the model.

```{r ch1_app_2_ppml_1}
tp_summary_app_2(
  formula = trade ~ 0 + log_dist_1986 + log_dist_1990 + log_dist_1994 +
    log_dist_1998 + log_dist_2002 + log_dist_2006 + cntg + lang + clny |
    exp_year + imp_year,
  data = filter(ch1_application2, importer != exporter),
  method = "ppml"
)
```

### Internal distance solution for the "distance puzzle"

This model requires us to add the internal distance variable to the PPML model and not filter the rows where the exporter and the importer are the same.

```{r ch1_app_2_intra_1}
tp_summary_app_2(
  formula = trade ~ 0 + log_dist_1986 + log_dist_1990 + log_dist_1994 +
    log_dist_1998 + log_dist_2002 + log_dist_2006 + cntg + lang + clny +
    log_dist_intra | exp_year + imp_year,
  data = ch1_application2,
  method = "ppml"
)
```

### Internal distance and home bias solution for the "distance puzzle"

This model requires us to add the same country variable to the internal distance model and repeat the rest of the steps from the last section.

```{r ch1_app_2_home_1}
tp_summary_app_2(
  formula = trade ~ log_dist_1986 + log_dist_1990 + log_dist_1994 +
    log_dist_1998 + log_dist_2002 + log_dist_2006 + cntg + lang + clny +
    log_dist_intra + smctry | exp_year + imp_year,
  data = ch1_application2,
  method = "ppml"
)
```

### Fixed effects solution for the "distance puzzle"

This model requires us to remove the internal distance and same country variables from the last model and include the internal pair variable to account for the intra-national fixed effects.

```{r ch1_app_2_fe_1}
tp_summary_app_2(
  formula = trade ~ 0 + log_dist_1986 + log_dist_1990 + log_dist_1994 +
    log_dist_1998 + log_dist_2002 + log_dist_2006 + cntg + lang + clny +
    intra_pair | exp_year + imp_year,
  data = ch1_application2,
  method = "ppml"
)
```

## Regional trade agreements effects

### Preparing the data

This model specification includes gravity covariates, including importer-time and exporter-time fixed effects, as in the equation

$$
\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 LANG_{i,j} +\right.\\
\text{ }& \:\left.\beta_4 CLNY_{i,j} + \beta_5 RTA_{ij,t}\right] \times \varepsilon_{ij,t}.
\end{align}
$$

In comparison to the previous examples, we need to create additional variables to include fixed effects that account for the observations where the exporter and the importer are the same. These variables are internal border, internal dyad and internal borders for different years.

The direct way of obtaining the desired variables is similar to what we did in the previous sections.

```{r ch1_app_3_data_1}
ch1_application3 <- agtpa_applications %>%
  filter(year %in% seq(1986, 2006, 4)) %>%
  mutate(
    exp_year = paste0(exporter, year),
    imp_year = paste0(importer, year),
    year = paste0("intl_border_", year),
    log_trade = log(trade),
    log_dist = log(dist),
    intl_brdr = ifelse(exporter == importer, pair_id, "inter"),
    intl_brdr_2 = ifelse(exporter == importer, 0, 1),
    pair_id_2 = ifelse(exporter == importer, "0-intra", pair_id)
  ) %>%
  pivot_wider(names_from = year, values_from = intl_brdr_2, values_fill = 0)
```

Notice that we used "0-intra" and not just "intra" because the rest of the observations in the internal dyads are numbers 1, ..., N, and R internals shall consider "0-intra" as the reference factor for being the first item when it orders the unique observations alphabetically. Also, observe the order of the resulting table, the pivoting of the table will put "0-intra" as the first row for the first exporter-importer dyad. This makes the difference between the expected or other behaviours in the next chapter.

In addition, we need to create the variable containing the trade sum to filter the cases where the sum by dyad is zero.

```{r ch1_app_3_data_2}
ch1_application3 <- ch1_application3 %>%
  group_by(pair_id) %>%
  mutate(sum_trade = sum(trade)) %>%
  ungroup()
```

### OLS standard RTA estimates with international trade only

The gravity specification, which includes $\pi_{i,t} + \chi_{i,t}$, means that we need to do something very similar to what we did in the last section.

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

```{r ch1_app_3_ols_1}
tp_summary_app_3(
  formula = log_trade ~ log_dist + cntg + lang + clny + rta | exp_year +
    imp_year,
  data = filter(ch1_application3, trade > 0, importer != exporter),
  method = "ols"
)
```

### PPML standard RTA estimates with international trade only

The model specification is very similar to OLS, and we only need to change the method specified in the function.

```{r ch1_app_3_ppml_1}
tp_summary_app_3(
  formula = trade ~ log_dist + cntg + lang + clny + rta | exp_year + imp_year,
  data = filter(ch1_application3, importer != exporter),
  method = "ppml"
)
```

### Addressing potential domestic trade diversion

The model specification is quite the same as PPML. We only need to add the international border variable but use the entire dataset instead of removing rows where the importer and the exporter are the same.

```{r ch1_app_3_intra_1}
tp_summary_app_3(
  formula = trade ~ log_dist + cntg + lang + clny + rta | exp_year + imp_year +
    intl_brdr,
  data = ch1_application3,
  method = "ppml"
)
```

### Addressing potential endogeneity of RTAs

The model specification includes the RTA variable and the exporter-time, importer-time and internal dyad fixed effects to account for domestic trade.

```{r ch1_app_3_endg_1}
tp_summary_app_3(
  formula = trade ~ rta | exp_year + imp_year + pair_id_2,
  data = filter(ch1_application3, sum_trade > 0),
  method = "ppml"
)
```

### Testing for potential "reverse causality" between trade and RTAs

We need to modify the previous model to include the forward lagged RTA variable (by four years) and consider where the trade sum is larger than zero.

```{r ch1_app_3_lead_1}
tp_summary_app_3(
  formula = trade ~ rta + rta_lead4 | exp_year + imp_year + pair_id_2,
  data = filter(ch1_application3, sum_trade > 0),
  method = "ppml"
)
```

### Addressing potential non-linear and phasing-in effects of RTAs

Instead of future-lagged RTA variables, as in the previous model, we modify the previous model and include the RTA backwards lagged variables instead.

```{r ch1_app_3_phsng_1}
tp_summary_app_3(
  formula = trade ~ rta + rta_lag4 + rta_lag8 + rta_lag12 | exp_year +
    imp_year + pair_id_2,
  data = filter(ch1_application3, sum_trade > 0),
  method = "ppml"
)
```

### Addressing globalization effects

In addition to the previous model, we include the international borders on different years besides the lagged RTAs.

```{r ch1_app_3_glbzn_1}
tp_summary_app_3(
  formula = trade ~ rta + rta_lag4 + rta_lag8 + rta_lag12 + intl_border_1986 +
    intl_border_1990 + intl_border_1994 + intl_border_1998 + intl_border_2002 |
    exp_year + imp_year + pair_id_2,
  data = filter(ch1_application3, sum_trade > 0),
  method = "ppml"
)
```

## References