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.
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:
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.
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, 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-16The package returns lm or glm objects instead of summaries. Doing that allows to use our functions in conjunction with broom or other packages.
S. L. 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.
S. L. 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 Scott L. Baier and Bergstrand (2009) and S. L. 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-16fit3 <- 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-16Silva 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
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -217.03 -28.10 -27.59 -23.42 1857.77
##
## 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: 8In 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
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -217.03 -28.10 -27.59 -23.42 1857.77
##
## 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: 8The 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-01The 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-01In 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-16In 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