# Basic model

Gravity models in their traditional form are inspired by Newton law of gravitation:

$F_{ij}=G\frac{M_{i}M_{j}}{D^{2}_{ij}}.$

The force $$F$$ between two bodies $$i$$ and $$j$$ with $$i \neq j$$ is proportional to the masses $$M$$ of these bodies and inversely proportional to the square of their geographical distance $$D$$. $$G$$ is a constant and as such of no major concern.

The underlying idea of a traditional gravity model, shown for international trade, is equally simple:

$X_{ij}=G\frac{Y_{i}^{\beta_{1}}Y_{j}^{\beta_{2}}}{D_{ij}^{\beta_{3}}}.$

The trade flow $$X$$ is explained by $$Y_{i}$$ and $$Y_{j}$$ that are the masses of the exporting and importing country (e.g. the GDP) and $$D_{ij}$$ that is the distance between the countries.

Dummy variables such as common borders $$contig$$ or regional trade agreements $$rta$$ can be added to the model. Let $$t_{ij}$$ be the transaction cost defined as:

$t_{ij}= D_{ij} \exp(contig_{ij} + rta_{ij})$

So that the model with friction becomes:

$X_{ij}=G\frac{Y_{i}^{\beta_{1}}Y_{j}^{\beta_{2}}}{t_{ij}^{\beta_{3}}}.$

A logarithmic operator can be applied to form a log-linear model and use a standard estimation methods such as OLS:

$\log X_{ij}=\beta_{0}\log G +\beta_{1}\log Y_{i}+\beta_{2}\log Y_{j}+\beta_{3}\log D_{ij}+\beta_{4}contig_{ij}+\beta_{5}rta_{ij}$

Provided trade barriers exists, the econometric literature proposes the Multilateral Resistance model defined by the equations:

$X_{ij}=\frac{Y_{i}Y_{j}}{Y}\frac{t_{ij}^{1-\sigma}}{P_{j}^{1-\sigma}\Pi_{i}^{1-\sigma}}$ with $P_{i}^{1-\sigma}=\sum_{j}\frac{t_{ij}^{1-\sigma}}{\Pi_{j}^{1-\sigma}}\frac{Y_{j}}{Y}.$ and $\Pi_{j}^{1-\sigma}=\sum_{i}\frac{t_{ij}^{1-\sigma}}{P_{i}^{1-\sigma}}\frac{Y_{i}}{Y}$

Basically the model proposes that the exports $$X_{ij}$$ from $$i$$ to $$j$$ are determined by the supply factors in $$i$$, $$Y_{i}$$, and the demand factors in $$j$$, $$Y_{j}$$, as well as the transaction costs $$t_{ij}$$.

Next to information on bilateral partners $$i$$ and $$j$$, information on the rest of the world is included in the gravity model with $$Y=\sum_{i} Y_{i}= \sum_{j} Y_{j}$$ that represents the worldwide sum of incomes (e.g. the world’s GDP).

In this model $$\sigma$$ represents the elasticity of substitution between all goods. A key assumption is to take a fixed value $$\sigma > 1$$ in order to account for the preference for a variation of goods (e.g. in this model goods can be replaced for other similar goods).

The Multilateral Resistance terms are included via the terms $$P$$, Inward Multilateral Resistance, and $$\Pi$$, Outward Multilateral Resistance. The Inward Multilateral Resistance $$P_i$$ is a function of the transaction costs of $$i$$ to all trade partners $$j$$. The Outward Multilateral Resistance $$\Pi_{j}$$ is a function of the transaction costs of $$j$$ to all trade partners $$i$$ and their demand.

The Multilateral Resistance terms dependent on each other. Hence, the estimation of structural gravity models becomes complex.

# Model estimation

To estimate gravity equations you need a square dataset including bilateral flows defined by the argument dependent_variable, a distance measure defined by the argument distance that is the key regressor, and other potential influences (e.g. contiguity and common currency) given as a vector in additional_regressors are required.

Some estimation methods require ISO codes or similar of type character variables to compute particular country effects. Make sure the origin and destination codes are of type “character”.

The rule of thumb for regressors or independent variables consists in:

• All dummy variables should be of type numeric (0/1).
• If an independent variable is defined as a ratio, it should be logged.

The user should perform some data cleaning beforehand to remove observations that contain entries that can distort estimates, notwithstanding the functions provided within this package will remove zero flows and distances.

See for gravity datasets and Stata code for estimating gravity models.

# Examples

All the examples here are adapted from Wölwer, Breßlein, and Burgard (2018). We included some examples that require further explanation as they perform some data transforming and therefore the functions provide a simplification for the end user.

## Double Demeaning

Double Demeaning, as introduced by Head and Mayer (2014), subtracts importer and exporter averages on the left and right hand side of the respective gravity equation, and all unilateral influences including the Multilateral Resistance terms vanish. Therefore, no unilateral variables may be added as independent variables for the estimation.

Our ddm function first logs the dependent variable and the distance variable. Afterwards, the dependent and independent variables are transformed in the following way (exemplary shown for trade flows, $$X_{ij}$$):

$(\log X_{ij})_{\text{DDM}} = (\log X_{ij}) - (\log X_{ij})_{\text{Origin Mean}} \\- (\log X_{ij})_{\text{Destination Mean}} + (\log X_{ij})_{\text{Mean}}.$

One subtracts the mean value for the origin country and the mean value for the destination country and adds the overall mean value to the logged trade flows. This procedure is repeated for all dependent and independent variables. The transformed variables are then used for the estimation.

DDM is easily applied, but, as shown in , its very sensitive to missing data.

An example of how to apply the function ddm to an example dataset in gravity and the resulting output is shown in the following:

library(gravity)

fit <- ddm(
dependent_variable = "flow",
distance = "distw",
code_origin = "iso_o",
code_destination = "iso_d",
data = gravity_no_zeros
)

summary(fit)
##
## Call:
## y_log_ddm ~ dist_log_ddm + rta_ddm + comcur_ddm + contig_ddm +
##     0
##
## Residuals:
##      Min       1Q   Median       3Q      Max
## -20.9411  -1.2268   0.3032   1.5195   8.4538
##
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## dist_log_ddm -1.60334    0.03312 -48.407   <2e-16 ***
## rta_ddm       0.79727    0.07004  11.383   <2e-16 ***
## comcur_ddm    0.17376    0.14613   1.189    0.234
## contig_ddm    1.00184    0.11981   8.362   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.318 on 17084 degrees of freedom
## Multiple R-squared:  0.2541, Adjusted R-squared:  0.254
## F-statistic:  1455 on 4 and 17084 DF,  p-value: < 2.2e-16

The package returns lm or glm objects instead of summaries. Doing that allows to use our functions in conjunction with broom or other packages, for example:

library(broom)

tidy(fit)
## # A tibble: 4 x 5
##   term         estimate std.error statistic  p.value
##   <chr>           <dbl>     <dbl>     <dbl>    <dbl>
## 1 dist_log_ddm   -1.60     0.0331    -48.4  0.
## 2 rta_ddm         0.797    0.0700     11.4  6.54e-30
## 3 comcur_ddm      0.174    0.146       1.19 2.34e- 1
## 4 contig_ddm      1.00     0.120       8.36 6.62e-17
glance(fit)
## # A tibble: 1 x 11
##   r.squared adj.r.squared sigma statistic p.value    df  logLik    AIC
##       <dbl>         <dbl> <dbl>     <dbl>   <dbl> <int>   <dbl>  <dbl>
## 1     0.254         0.254  2.32     1455.       0     4 -38611. 77232.
## # … with 3 more variables: BIC <dbl>, deviance <dbl>, df.residual <int>
augment(fit)
## # A tibble: 17,088 x 12
##    y_log_ddm dist_log_ddm rta_ddm comcur_ddm contig_ddm  .fitted .se.fit
##        <dbl>        <dbl>   <dbl>      <dbl>      <dbl>    <dbl>   <dbl>
##  1  -0.710         0.540   0.0175    0.0562    -0.0435  -0.885   0.0205
##  2   0.00279      -0.0201  0.0318    0.0501    -0.00474  0.0616  0.00754
##  3  -0.373         0.154  -0.197    -0.00870   -0.0502  -0.456   0.0126
##  4  -3.85         -0.758  -0.0408    0.0562    -0.0533   1.14    0.0298
##  5   0.301         0.223  -0.191    -0.00705   -0.0237  -0.534   0.0111
##  6  -3.45          0.599  -0.0747   -0.0747    -0.0762  -1.11    0.0176
##  7   1.80         -0.861   0.0562    0.0562    -0.0152   1.42    0.0276
##  8  -1.43         -0.0245 -0.127     0.0562    -0.0319  -0.0845  0.0135
##  9   3.10         -0.0260 -0.0195    0.0562    -0.0426  -0.00683 0.0104
## 10  -3.88          0.622  -0.0227    0.0562    -0.0486  -1.05    0.0217
## # … with 17,078 more rows, and 5 more variables: .resid <dbl>, .hat <dbl>,
## #   .sigma <dbl>, .cooksd <dbl>, .std.resid <dbl>

## Bonus Vetus

Baier and Bergstrand (2010) suggests a modification of the simple OLS that is easily implemented, allows for comparative statics and yields results close to those of NLS, called Bonus vetus OLS (BVU and BVW). They estimate gravity models in their additive form.

As unilateral income elasticities are assumed, flows are divided by the product of unilateral incomes. The dependent variable for the estimation is therefore $\log\left(\frac{y}{inc_{o} \: inc_{d}}\right).$

By applying a Taylor-series expansion and the assumption of symmetrical, bilateral trade costs, they develop a reduced gravity model in which multilateral and worldwide resistance enter exogenously.

Baier and Bergstrand (2010) distinguishes two types of Bonus vetus estimations depending on how the Taylor-series is centered. One method, called BVU, uses simple averages while the other, called BVW, uses GDP weights. Depending on which method is used, the transaction costs are weighted differently. For advantages and disadvantages of both methods see Baier and Bergstrand (2009) and Baier and Bergstrand (2010).

To give an example with simple averages (BVU), distance would be transformed to Multilateral and World Resistance in the following way: $MWR_{ij} = \frac{1}{N}\left(\sum_{i=1}^{N}\log D_{ij} \right)+\frac{1}{N}\left(\sum_{j=1}^{N}\log D_{ij} \right)-\frac{1}{N^{2}}\left(\sum_{i=1}^{N}\sum_{j=1}^{N}\log D_{ij} \right)$ with $$D_{ij}$$ denoting the bilateral distance, $$N$$ the number of countries and $$(MWR)_{D_{ij}}$$ the transformed variable adjusted for multilateral resistances.

When using weighted averages (BVW), the simple averages are replaced by GDP weights. The transformed variables are included as independent variables in the estimation. The resulting equation can be estimated by simple OLS.

An example of how to apply the functions bvu and bvw to an example dataset in gravity and the resulting output is shown in the following:

fit2 <- bvu(
dependent_variable = "flow",
distance = "distw",
income_origin = "gdp_o",
income_destination = "gdp_d",
code_origin = "iso_o",
code_destination = "iso_d",
data = gravity_no_zeros
)

tidy(fit2)
## # A tibble: 5 x 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)  -20.2      0.0192  -1051.   0.
## 2 dist_log_mr   -1.69     0.0359    -47.0  0.
## 3 rta_mr         0.582    0.0793      7.35 2.12e-13
## 4 contig_mr      1.00     0.131       7.62 2.75e-14
## 5 comcur_mr      0.191    0.167       1.14 2.53e- 1
fit3 <- bvw(
dependent_variable = "flow",
distance = "distw",
income_origin = "gdp_o",
income_destination = "gdp_d",
code_origin = "iso_o",
code_destination = "iso_d",
data = gravity_no_zeros
)

tidy(fit3)
## # A tibble: 5 x 5
##   term        estimate std.error statistic   p.value
##   <chr>          <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)  -20.4      0.0546   -373.   0.
## 2 dist_log_mr   -0.634    0.0257    -24.7  1.09e-132
## 3 rta_mr         1.31     0.0670     19.5  4.07e- 84
## 4 comcur_mr     -0.896    0.142      -6.31 2.86e- 10
## 5 contig_mr      1.38     0.126      10.9  1.11e- 27

## Poisson Pseudo Maximum Likelihood (PPML)

Silva and Tenreyro (2006) argue that estimating gravity equations in their additive form by OLS leads to inconsistency in the presence of heteroscedasticity and advice to estimate gravity models in their multiplicative form.

An example of how to apply the function ppml to an example dataset in gravity and the resulting output is shown in the following:

fit4 <- ppml(
dependent_variable = "flow",
distance = "distw",
data = gravity_no_zeros
)

tidy(fit4)
## # A tibble: 5 x 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)   5.76      0.772      7.47  8.54e-14
## 2 dist_log      0.0243    0.0873     0.278 7.81e- 1
## 3 rta           1.27      0.170      7.46  9.20e-14
## 4 comcur        1.11      0.180      6.15  7.89e-10
## 5 contig        1.76      0.181      9.69  3.71e-22

In order to obtain robust standard errors (i.e. in a similar way to vce(robust) in Stata) you can include robust = T to the arguments:

fit4 <- ppml(
dependent_variable = "flow",
distance = "distw",
robust = TRUE,
data = gravity_no_zeros
)

tidy(fit4)
## # A tibble: 5 x 4
##   estimate std.error statistic  p.value
##      <dbl>     <dbl>     <dbl>    <dbl>
## 1   5.76      0.772      7.47  8.54e-14
## 2   0.0243    0.0873     0.278 7.81e- 1
## 3   1.27      0.170      7.46  9.20e-14
## 4   1.11      0.180      6.15  7.89e-10
## 5   1.76      0.181      9.69  3.71e-22

## Gamma Pseudo Maximum Likelihood (GPML)

The estimation method is similar to PPML, but utilizes the gamma instead of the poisson distribution, thereby implies different assumptions to the data structure and does not allow for zero trade values.

Silva and Tenreyro (2006) argue in favor of PPML instead of GPML, especially in case of heteroscedasticity, Head and Mayer (2014) highlight that depending on data structure there exist cases where GPML is preferable to PPML.

An example of how to apply the function gpml to an example dataset in gravity and the resulting output is shown in the following:

fit5 <- gpml(
dependent_variable = "flow",
distance = "distw",
robust = TRUE,
data = gravity_no_zeros
)

tidy(fit5)
## # A tibble: 5 x 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)    7.40      1.03       7.20 5.81e-13
## 2 dist_log      -0.161     0.118     -1.36 1.73e- 1
## 3 rta            1.03      0.160      6.42 1.32e-10
## 4 comcur         0.783     0.183      4.28 1.88e- 5
## 5 contig         1.71      0.297      5.74 9.49e- 9

## Negative Binomial Pseudo Maximum Likelihood (NBPML)

The estimation method is similar to PPML, but utilizes the negative binomial instead of the poisson distribution, thereby implies different assumptions to the data structure and does not allow for zero trade values.

An example of how to apply the function nbpml to an example dataset in gravity and the resulting output is shown in the following:

fit6 <- nbpml(
dependent_variable = "flow",
distance = "distw",
robust = TRUE,
data = gravity_no_zeros
)

tidy(fit6)
## # A tibble: 5 x 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)    7.40      1.03       7.20 5.81e-13
## 2 dist_log      -0.161     0.118     -1.36 1.73e- 1
## 3 rta            1.03      0.160      6.43 1.32e-10
## 4 comcur         0.783     0.183      4.28 1.89e- 5
## 5 contig         1.71      0.297      5.74 9.51e- 9

In order to use the fixed effects method with panel data, a huge number of dummy variables has to be included into the estimation. Thus, estimating these models can lead to significant computational difficulties.

Head, Mayer, and Ries (2010) present Tetrads as an estimation method circumventing this problem. They exploit the multiplicative form of the gravity equation to form the ratio of ratios. In doing so, both MR terms drop out of the equation.

An example of how to apply the function tetrads to an example dataset in gravity and the resulting output is shown in the following:

fit8 <- tetrads(
dependent_variable = "flow",
distance = "distw",
code_origin = "iso_o",
code_destination = "iso_d",
filter_origin = "CHN",
filter_destination = "USA",
data = gravity_no_zeros
)

tidy(fit8)
## # A tibble: 3 x 5
##   term         estimate std.error statistic  p.value
##   <chr>           <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)    -0.400    0.0216     -18.5 1.47e-75
## 2 dist_log_rat   -1.32     0.0231     -57.4 0.
## 3 rta_rat         0.883    0.0557      15.9 2.80e-56

In addition to robust standard errors as in the previous examples, in the case of tetrads you can also be interested in computing multiway variance-covariance, an it can be done by adding multiway = T to the arguments:

fit8 <- tetrads(
dependent_variable = "flow",
distance = "distw",
code_origin = "iso_o",
code_destination = "iso_d",
filter_origin = "CHN",
filter_destination = "USA",
multiway = TRUE,
data = gravity_no_zeros
)

tidy(fit8)
## # A tibble: 3 x 5
##   term         estimate std.error statistic  p.value
##   <chr>           <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)    -0.400    0.0215     -18.6 3.05e-76
## 2 dist_log_rat   -1.32     0.0247     -53.7 0.
## 3 rta_rat         0.883    0.0487      18.1 1.11e-72