A crash course on gravity models

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}$

Trade barriers model

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",
  additional_regressors = c("rta", "comcur", "contig"),
  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.

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",
  additional_regressors = c("rta", "contig", "comcur"),
  income_origin = "gdp_o",
  income_destination = "gdp_d",
  code_origin = "iso_o",
  code_destination = "iso_d",
  data = gravity_no_zeros
)

summary(fit2)

## 
## Call:
## y_log_bvu ~ dist_log_mr + rta_mr + contig_mr + comcur_mr
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -20.0642  -1.2305   0.2932   1.5810   9.6659 
## 
## Coefficients:
##              Estimate Std. Error   t value Pr(>|t|)    
## (Intercept) -20.15771    0.01919 -1050.573  < 2e-16 ***
## dist_log_mr  -1.68697    0.03586   -47.048  < 2e-16 ***
## rta_mr        0.58230    0.07926     7.347 2.12e-13 ***
## contig_mr     1.00135    0.13148     7.616 2.75e-14 ***
## comcur_mr     0.19102    0.16727     1.142    0.253    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.508 on 17083 degrees of freedom
## Multiple R-squared:  0.2294, Adjusted R-squared:  0.2292 
## F-statistic:  1271 on 4 and 17083 DF,  p-value: < 2.2e-16

fit3 <- bvw(
  dependent_variable = "flow",
  distance = "distw",
  additional_regressors = c("rta", "comcur", "contig"),
  income_origin = "gdp_o",
  income_destination = "gdp_d",
  code_origin = "iso_o",
  code_destination = "iso_d",
  data = gravity_no_zeros
)

summary(fit3)

## 
## Call:
## y_log_bvw ~ dist_log_mr + rta_mr + comcur_mr + contig_mr
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -19.1204  -1.3366   0.4068   1.7068  10.6927 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -20.37260    0.05458 -373.27  < 2e-16 ***
## dist_log_mr  -0.63445    0.02566  -24.73  < 2e-16 ***
## rta_mr        1.30989    0.06703   19.54  < 2e-16 ***
## comcur_mr    -0.89609    0.14201   -6.31 2.86e-10 ***
## contig_mr     1.37696    0.12606   10.92  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.683 on 17083 degrees of freedom
## Multiple R-squared:  0.1185, Adjusted R-squared:  0.1183 
## F-statistic:   574 on 4 and 17083 DF,  p-value: < 2.2e-16

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",
  additional_regressors = c("rta", "comcur", "contig"),
  data = gravity_no_zeros
)

summary(fit4)

## 
## Call:
## y_ppml ~ dist_log + rta + comcur + contig
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.76232    0.77160   7.468 8.54e-14 ***
## dist_log     0.02428    0.08728   0.278    0.781    
## rta          1.26582    0.16972   7.458 9.20e-14 ***
## comcur       1.10604    0.17982   6.151 7.89e-10 ***
## contig       1.75821    0.18141   9.692  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 42366.31)
## 
##     Null deviance: 78081855  on 17087  degrees of freedom
## Residual deviance: 61989121  on 17083  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 8

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:

fit4r <- ppml(
  dependent_variable = "flow",
  distance = "distw",
  additional_regressors = c("rta", "comcur", "contig"),
  robust = TRUE,
  data = gravity_no_zeros
)

summary(fit4r)

## 
## Call:
## y_ppml ~ dist_log + rta + comcur + contig
## 
## Coefficients:
##      Estimate Std. Error t value Pr(>|t|)    
## [1,]  5.76232    0.77160   7.468 8.54e-14 ***
## [2,]  0.02428    0.08728   0.278    0.781    
## [3,]  1.26582    0.16972   7.458 9.20e-14 ***
## [4,]  1.10604    0.17982   6.151 7.89e-10 ***
## [5,]  1.75821    0.18141   9.692  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 42366.31)
## 
##     Null deviance: 78081855  on 17087  degrees of freedom
## Residual deviance: 61989121  on 17083  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 8

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",
  additional_regressors = c("rta", "comcur", "contig"),
  robust = TRUE,
  data = gravity_no_zeros
)

summary(fit5)

##     Estimate         Std. Error        z value          Pr(>|z|)        
##  Min.   :-0.1613   Min.   :0.1183   Min.   :-1.364   Min.   :0.000e+00  
##  1st Qu.: 0.7833   1st Qu.:0.1603   1st Qu.: 4.279   1st Qu.:0.000e+00  
##  Median : 1.0295   Median :0.1831   Median : 5.740   Median :1.000e-08  
##  Mean   : 2.1508   Mean   :0.3571   Mean   : 4.457   Mean   :3.453e-02  
##  3rd Qu.: 1.7065   3rd Qu.:0.2973   3rd Qu.: 6.424   3rd Qu.:1.879e-05  
##  Max.   : 7.3963   Max.   :1.0266   Max.   : 7.205   Max.   :1.726e-01

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",
  additional_regressors = c("rta", "comcur", "contig"),
  robust = TRUE,
  data = gravity_no_zeros
)

summary(fit6)

##     Estimate         Std. Error        z value          Pr(>|z|)        
##  Min.   :-0.1613   Min.   :0.1183   Min.   :-1.364   Min.   :0.000e+00  
##  1st Qu.: 0.7831   1st Qu.:0.1603   1st Qu.: 4.277   1st Qu.:0.000e+00  
##  Median : 1.0297   Median :0.1831   Median : 5.739   Median :1.000e-08  
##  Mean   : 2.1508   Mean   :0.3571   Mean   : 4.457   Mean   :3.454e-02  
##  3rd Qu.: 1.7064   3rd Qu.:0.2973   3rd Qu.: 6.425   3rd Qu.:1.893e-05  
##  Max.   : 7.3961   Max.   :1.0266   Max.   : 7.205   Max.   :1.727e-01

Tetrads

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",
  additional_regressors = "rta",
  code_origin = "iso_o",
  code_destination = "iso_d",
  filter_origin = "CHN",
  filter_destination = "USA",
  data = gravity_no_zeros
)

summary(fit8)

## 
## Call:
## y_log_rat ~ dist_log_rat + rta_rat
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -21.3989  -1.2542   0.3247   1.5871  10.3864 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -0.39969    0.02162  -18.49   <2e-16 ***
## dist_log_rat -1.32417    0.02307  -57.40   <2e-16 ***
## rta_rat       0.88305    0.05566   15.87   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.77 on 16710 degrees of freedom
## Multiple R-squared:  0.2537, Adjusted R-squared:  0.2536 
## F-statistic:  2840 on 2 and 16710 DF,  p-value: < 2.2e-16

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",
  additional_regressors = "rta",
  code_origin = "iso_o",
  code_destination = "iso_d",
  filter_origin = "CHN",
  filter_destination = "USA",
  multiway = TRUE,
  data = gravity_no_zeros
)

summary(fit8)

##     Estimate         Std. Error         t value            Pr(>|t|)        
##  Min.   :-1.3242   Min.   :0.02152   Min.   :-53.6833   Min.   :0.000e+00  
##  1st Qu.:-0.8619   1st Qu.:0.02309   1st Qu.:-36.1287   1st Qu.:1.525e-76  
##  Median :-0.3997   Median :0.02467   Median :-18.5741   Median :3.049e-76  
##  Mean   :-0.2803   Mean   :0.03164   Mean   :-18.0460   Mean   :3.701e-73  
##  3rd Qu.: 0.2417   3rd Qu.:0.03670   3rd Qu.: -0.2273   3rd Qu.:5.552e-73  
##  Max.   : 0.8830   Max.   :0.04873   Max.   : 18.1195   Max.   :1.110e-72

References

Baier, Scott L., and Jeffrey H. Bergstrand. 2009. “Bonus Vetus Ols: A Simple Method for Approximating International Trade-Cost Effects Using the Gravity Equation.” Journal of International Economics 77 (1): 77–85. https://doi.org/10.1016/j.jinteco.2008.10.004.

Baier, S. L., and J. H. Bergstrand. 2010. “The Gravity Model in International Trade: Advances and Applications.” In, edited by Peter A. G. van Bergeijk and StevenEditors Brakman. Cambridge University Press. https://doi.org/10.1017/CBO9780511762109.

Head, Keith, and Thierry Mayer. 2014. “Chapter 3 - Gravity Equations: Workhorse,Toolkit, and Cookbook.” In Handbook of International Economics, edited by Gita Gopinath, Elhanan Helpman, and Kenneth Rogoff, 4:131–95. Handbook of International Economics. Elsevier. https://doi.org/10.1016/B978-0-444-54314-1.00003-3.

Head, Keith, Thierry Mayer, and John Ries. 2010. “The Erosion of Colonial Trade Linkages After Independence.” Journal of International Economics 81 (1): 1–14. https://doi.org/10.1016/j.jinteco.2010.01.002.

Silva, J. M. C. Santos, and Silvana Tenreyro. 2006. “The Log of Gravity.” The Review of Economics and Statistics 88 (4): 641–58. https://doi.org/10.1162/rest.88.4.641.

Wölwer, Anna-Lena, Martin Breßlein, and Jan Pablo Burgard. 2018. “Gravity Models in R.” Austrian Journal of Statistics 47 (4): 16–35. https://doi.org/10.17713/ajs.v47i4.688.

Mauricio Vargas

2024-11-21