伝統的な重力方程式の推定
田中 鮎夢
2023-12-01
This version: 2025-12-16
1 貿易データ
本ページでは、重力方程式の基本的な推定方法を説明する。日本語の文献としては、重力方程式について、伊藤・田中『現実からまなぶ国際経済学』 や 田中「国際貿易と重力の意外な関係: 重力方程式の基本」 に説明がある。
gravity.xlsx データには、2005 年の 1年間の世界の2国間の貿易額が収録されている。
Rでデータを読み込む。
library(readxl)
gravity <- read_excel("gravity.xlsx")
head(gravity)## # A tibble: 6 × 20
## importer exporter year imports gdp_exporter gdp_importer join_exporter
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 AFG ABW 2005 0 NA 6814753581 NA
## 2 AGO ABW 2005 0 NA 30632364954 NA
## 3 ALB ABW 2005 0 NA 8376483740 NA
## 4 ANT ABW 2005 4335. NA NA NA
## 5 ARE ABW 2005 0.951 NA 133000000000 NA
## 6 ARG ABW 2005 0.658 NA 183000000000 NA
## # ℹ 13 more variables: join_importer <dbl>, exporternum <dbl>,
## # importernum <dbl>, contig <dbl>, comlang_off <dbl>, colony <dbl>,
## # dist <dbl>, REPlandlocked <dbl>, PARTlandlocked <dbl>, religion <dbl>,
## # onein <dbl>, bothin <dbl>, nonein <dbl>オンラインから直接データを読み込みたい場合は、パッケージgdataをインストールした上で、以下のコードを用いる。
# install.packages("gdata")
library(gdata)
library(readxl)
download.file("https://ayumu-tanaka.github.io/teaching/gravity.xlsx","gravity.xlsx",mode="wb")
gravity <- read_excel("gravity.xlsx")変数 importer は輸入国の ISO コードと呼ばれる 3桁の国コードである。同様に、exporter は輸出国 の ISO コードである。例えば、日本の ISO コードは JPN、アメリカの ISO コードは USA である。また、ISOコード（exporter, importer）の他に、輸入額（imports）、輸出国のGDP（gdp_exporter）、輸入国のGDP（gdp_importer）、輸出国と輸入国の間の距離（dist）が含まれている。
2 対数変換
R で次のような重力方程式を推定することを考える。
\[ 輸入額 = 定数 \times \frac{輸出国のGDP^{\alpha} \times 輸出国のGDP^{\beta}}{輸出国と輸入国の間の距離^{\gamma}} \times e ^{\delta \times 言語の共通性} \]
この式は非線形のためコンピューターの能力上、推定が難しかった。そのため、従来は両辺の対数をとって、線形にされてきた。線形にされた式の推定は最小二乗法で比較的簡単に推定できる。そこでまず、回帰分析に用いる変数の対数を取る。Rのコードは次の通りである。
gravity$limports <- log(gravity$imports)
gravity$lgdp_exporter <- log(gravity$gdp_exporter)
gravity$lgdp_importer <- log(gravity$gdp_importer)
gravity$ldist <- log(gravity$dist)ここで、0より大きい値しか対数値にできないことに注意が必要である。 そのため、輸入額が0より大きい値のサンプルを作成しておく。
gravity2 <- subset(gravity,gravity$imports>0)3 推定
通常、貿易データは不均一分散の性質を持つので、デフォルトのlmで推定するのではなく、feolsで推定する。そのために、パッケージfixestをインストールする。
install.packages("fixest")パッケージfixestの使い方は、公式ページで確認できる。
不均一分散頑健な標準誤差を計算するよう、feolsで重力方程式を推定するコードは以下の通りである。summary(ols)により回帰分析の結果が出力される。
library(fixest)
ols <- fixest::feols(limports ~ lgdp_exporter + lgdp_importer
+ ldist + comlang_off,
data = gravity2,
vcov = "hetero")
fixest::etable(ols)## ols
## Dependent Var.: limports
##
## Constant -33.75*** (0.3529)
## lgdp_exporter 1.226*** (0.0076)
## lgdp_importer 0.9510*** (0.0076)
## ldist -1.374*** (0.0204)
## comlang_off 1.293*** (0.0504)
## _______________ __________________
## S.E. type Heteroskedas.-rob.
## Observations 19,978
## R2 0.64238
## Adj. R2 0.64230
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14 推定結果の解釈
4.1 連続変数
推定結果において、たとえば、輸出国のGDP (lgdp_exporter) の推定係数は\(\hat{\alpha}=1.226\)になっている。これは、以下のような関係にあると言える。
\[ \frac{\partial \hat{\mathrm{ln}輸入額}}{\partial \mathrm{ln}輸出国のGDP} = \hat{\alpha}=1.226 \]
また、対数微分の公式 \[ \frac{\mathrm{dln}x}{\mathrm{d}x}=\frac{1}{x} \] から、\(\mathrm{dln}x=\frac{\mathrm{d}x}{x}\)となるので、\(\mathrm{dln}x\)は変化率を表している。この関係を使って上の式を書き直すと、近似的に以下が成り立つ。
\[ \frac{輸入額の変化率}{輸出国のGDPの変化率} = \hat{\alpha}=1.226 \]
この式は、弾性値（=Yの変化率/Xの変化率）の形になっている。そのため、輸出国のGDPが1%大きくなれば、平均的に輸入額が1.226%大きくなる傾向にあると表現できる。
4.2 ダミー変数
また、言語の共通性 (comlang_off) の推定係数は\(\hat{\delta}=1.293\)である。もともと、言語の共通性以外の他の項を\(B\)としてまとめると、重力方程式の定式化から以下のような関係にある。
\[ 輸入額 = B \times e ^{\delta \times 言語の共通性} \] そのため、言語が共通の場合は、 \[ \hat{輸入額} = B \times e ^{1.293 \times 1} = B \times e ^{1.293} \] そのため、言語が共通ではない場合は、
\[ \hat{輸入額} = B \times e ^{1.293 \times 0} = B \times e ^{0} = B \]
となる。言語が共通な場合は、
\[ \frac{B \times e ^{1.293}}{B}=e ^{1.293} \fallingdotseq 3.64 \]
より、平均的に3.64倍貿易額が大きい傾向にあることが分かる。
なお、以下のように計算している。
b <- exp(1.293)
b## [1] 3.6437015 推定結果の整形と表示1
with modelsummary
fixestパッケージのetable関数を用いて推定結果を表示できるが、より見栄えの良い推定結果表を作成したい場合もある。ここでは、Arel-Bundock (2022)
が開発したmodelsummaryを用いて、推定結果を表示する。
#install.packages("modelsummary")
library(modelsummary)
#(1) lmを使用して通常の回帰分析を実行し、ols1に保存。
ols1 <- lm(limports ~ lgdp_exporter
+ lgdp_importer + ldist + comlang_off,
data = gravity2)
#(2) lm_robustを使用して回帰分析を実行し、不均一分散頑健な標準誤差を計算し、ols2に保存。
ols2 <- fixest::feols(limports ~ lgdp_exporter
+ lgdp_importer + ldist + comlang_off,
vcov = "hetero",
data = gravity2)
modelsummary(list(ols1, ols2), stars = TRUE)| (1) | (2) | |
|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | ||
| (Intercept) | -33.751*** | -33.751*** |
| (0.350) | (0.353) | |
| lgdp_exporter | 1.226*** | 1.226*** |
| (0.008) | (0.008) | |
| lgdp_importer | 0.951*** | 0.951*** |
| (0.008) | (0.008) | |
| ldist | -1.374*** | -1.374*** |
| (0.022) | (0.020) | |
| comlang_off | 1.293*** | 1.293*** |
| (0.051) | (0.050) | |
| Num.Obs. | 19978 | 19978 |
| R2 | 0.642 | 0.642 |
| R2 Adj. | 0.642 | 0.642 |
| AIC | 93181.8 | 93179.8 |
| BIC | 93229.2 | 93219.3 |
| Log.Lik. | -46584.884 | |
| RMSE | 2.49 | 2.49 |
| Std.Errors | Heteroskedasticity-robust | |
より手の込んだ推定結果表も作成できる。 参考：modelsummary: regression tables
modelsummary(
list(ols1, ols2),
estimate = c("estimate",
"{estimate}{stars}"),
coef_omit = "Intercept",
gof_omit = 'DF|Deviance|AIC|BIC|RMSE',
notes = list('* p<0.1; ** p<0.05; *** p<0.01')
)| (1) | (2) | |
|---|---|---|
| * p<0.1; ** p<0.05; *** p<0.01 | ||
| lgdp_exporter | 1.226 | 1.226*** |
| (0.008) | (0.008) | |
| lgdp_importer | 0.951 | 0.951*** |
| (0.008) | (0.008) | |
| ldist | -1.374 | -1.374*** |
| (0.022) | (0.020) | |
| comlang_off | 1.293 | 1.293*** |
| (0.051) | (0.050) | |
| Num.Obs. | 19978 | 19978 |
| R2 | 0.642 | 0.642 |
| R2 Adj. | 0.642 | 0.642 |
| Log.Lik. | -46584.884 | |
| Std.Errors | Heteroskedasticity-robust | |
6 推定結果のLaTex形式出力
多くの経済学者が、Microsoft Wordではなく、LaTexを用いて論文を執筆している。最近では、LaTexをオンラインで用いることができるOverleaf (Online LaTeX Editor)が人気である。Overleafは、無料でアカウントを作成して利用できる。Posit Cloudで作成した推定結果表をOverleafにアップロードして、論文に組み込む作業を考えてみる。
まず、modelsummaryを用いて、以下のコードのように推定結果をLaTex形式で出力することができる。
modelsummary(list(ols1, ols2), output = "02_table.tex")生成された”table.tex”をOverleaf にアップロードして、LaTex形式で表を作成することができる。なお、LaTexのpreambleに以下のコード記載することが必要である。
\usepackage{tabularray}
\UseTblrLibrary{booktabs}
\usepackage{longtable}
\usepackage{array}
\usepackage{multirow}
\usepackage{wrapfig} % if necessary
\usepackage{float}
\usepackage{colortbl} % if necessary
\usepackage{pdflscape}
\usepackage{tabu}
\usepackage{threeparttable}
\usepackage{threeparttablex}
\usepackage[normalem]{ulem}
\usepackage{makecell}
\usepackage{xcolor} % if necessary
\usepackage{siunitx}上記ののpreambleは、modelsummaryのバージョン2.0.0の場合のものである。
Overleafのサンプルページはこちらにある。
なお、英語の文書の場合、RマークダウンをKnit to PDF により、PDFファイルに変換する過程で、LaTex形式に変換されている。そのため、RマークダウンをPDFに変換することで、LaTex形式の表を生成することもできる。しかし、日本語の文書の場合、Knit to PDFでは、文字コードに起因するエラーが発生して、対処が難しい。
Overleafを使うことで、Posit Cloudで生成したLaTex形式の表を含む日本語の論文を作成することも比較的容易にできる。
7 係数プロット with
modelplot
library(ggplot2)
modelplot(list(ols1, ols2), coef_omit = "Intercept")
ggsave("02_coefplot.png")補論：データの出所
UN and WTO. “A Practical Guide to Trade Policy Analysis,” Chapter 3 で用いられているデータ「gravity.dta」をExcelに変換したデータを本文書では用いた。元データは以下のサイトより無料で入手できる。
本文書で用いたデータのStata版「gravity.dta」は上記のサイトで提供されているdo-file「BuildingDatabaseApproach.do」を実行することで生成される。有斐閣の教科書サポートページでも、データ「gravity.dta」を提供しているので、ダウンロードして用いることもできる。
補論：estimatr
Posit Cloudでは、estimatrパッケージ（Blair et al. (2018)
が開発）のインストールに失敗する場合がある。そのため、estimatrパッケージではなく、fixestパッケージを用いる。
estimatrパッケージを用いた不均一分散頑健な重力方程式の推定コードは以下の通りである。
library(estimatr)
ols <- lm_robust(limports ~
lgdp_exporter + lgdp_importer + ldist + comlang_off,
data = gravity2)