Nothing
A crash course on gravity models
In gravity: Estimation Methods for Gravity Models
knitr::opts_chunk$set(eval = TRUE, message = FALSE, warning = FALSE)
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 \href{https://sites.google.com/site/hiegravity/}{Gravity Equations: Workhorse, Toolkit, and Cookbook}
for gravity datasets and Stata code for estimating gravity models.
Examples
All the examples here are adapted from @WoelwerBressleinBurgard2018. 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 @Head2014, 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 \cite{Head2014}, 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)
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
@Baier2010 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.
@Baier2010 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 @Baier2009 and @Baier2010.
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)
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)
Poisson Pseudo Maximum Likelihood (PPML)
@Santos2006 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)
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)
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.
@Santos2006 argue in favor of PPML instead of GPML, especially in case of heteroscedasticity, @Head2014 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)
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)
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.
@Head2010 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)
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)
References
\insertRef{Head2014}{gravity}
\insertRef{WoelwerBressleinBurgard2018}{gravity}
Try the gravity package in your browser
Any scripts or data that you put into this service are public.
gravity documentation built on May 2, 2023, 9:13 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(eval = TRUE, message = FALSE, warning = FALSE)
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 \href{https://sites.google.com/site/hiegravity/}{Gravity Equations: Workhorse, Toolkit, and Cookbook} for gravity datasets and Stata code for estimating gravity models.
Examples
All the examples here are adapted from @WoelwerBressleinBurgard2018. 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 @Head2014, 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 \cite{Head2014}, 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)
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
@Baier2010 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.
@Baier2010 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 @Baier2009 and @Baier2010.
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)
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)
Poisson Pseudo Maximum Likelihood (PPML)
@Santos2006 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)
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)
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.
@Santos2006 argue in favor of PPML instead of GPML, especially in case of heteroscedasticity, @Head2014 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)
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)
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.
@Head2010 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)
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)
References
\insertRef{Head2014}{gravity} \insertRef{WoelwerBressleinBurgard2018}{gravity}
Try the gravity package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.