This version: 2024-11-17

1 はじめに

本ページでは、政策が２国間の貿易額に与える影響を分析するため、2期間２グループの「2 x 2」の差の差 (DiD: difference-in-differences) 分析を実施する。単純なOLSによるDiD分析から始めて、 Anderson and Van Wincoop (2003) が指摘した多角的貿易抵抗指数（物価効果）を制御するために、時変の輸出国（輸入国）固定効果、貿易ペアの固定効果、さらには年次固定効果を加えて、重力方程式をポワソン擬似最尤（PPML）で推定していく。

本ページでは、1958年1月1日に、ローマ条約が発効し、欧州経済共同体（EEC、のちの欧州共同体 (EC)）が発足した事例を用いて、EEC発足が２国間貿易に与えた影響を推定する。EEC発足は、1958年という単一時点のため、本ページで扱う処置は時差のない処置である。時差のある処置の場合は、処置効果の推定はより難しくなるが、本ページで扱う時差のない処置は、そうした困難はない（Goodman-Bacon (2021), Roth et al. (2023)）。

推定には、本ページでは、Berge and McDermott (2023) によって開発されたfixestパッケージを用いる。パッケージのインストールは以下のコードで行われる。

install.packages("fixest")

パッケージの読み込みを行う。

# パッケージの読み込み
library(fixest)

2 貿易データ

gravityIMFeu.dtaには、フランスの研究機関CEPIIが作成している、Gravityという無料のデータベースから、２国間の貿易額やRTAの有無など必要最小限のデータを抜き出している。gravityIMFeu.dtaには、1948–1972年の25年間の世界の2国間の貿易額 (tradeflow_imf_d, tradeflow_imf_o) が収録されている。今回は古い期間の貿易データを用いるため、1996年からしかデータがないCEPIIが作成した貿易データ(tradeflow_baci) ではなく、IMFの貿易データ(tradeflow_imf_d) を用いる。また、一般に、関税計算のために目的地の貿易統計の方が正確であると言われているので、目的地の貿易額 (tradeflow_imf_d) を分析に用いる。

Gravityデータベースは定期的に更新されている。本ページで用いるデータは、「Gravity_dta_V202211」版であり、2024年3月に入手した。主な変数は以下の通りである。

year: Year
rta: 1 = RTA (source: WTO)
pair: group(iso3_o iso3_d)
iso3_o: Origin ISO3 alphabetic
iso3_d: Destination ISO3 alphabetic
gdp_o: Origin GDP (current thousands US$)
gdp_d: Destination GDP (current thousands US$)
gdpcap_o: Origin GDP per cap (current thousands US$)
gdpcap_d: Destination GDP per cap (current thousands US$)
treated_after: 1 = Both Origin and Destination EU members and after 1958
treated: 1 = Both Origin and Destination EU members
after: 1 = after 1958
tradeflow_imf_o: Trade flows as reported by the origin, 1000
current USD (source: IMF)
tradeflow_imf_d: Trade flows as reported by the destination, 1000 current USD (source: IMF)

RでStata形式のデータ（.dta）を読み込むに、パッケージrioを用いる。rioは、Stata形式に限らず、様々な拡張子のデータをimport()関数で読み込むことができる。以下でインストールする。

install.packages("rio")

そして、rioパッケージのimport関数でdtaファイルを読み込む。本来、library(rio)でrioパッケージを読み込めば、import関数の前にrio::をつける必要はないが、どのパッケージを用いているのかを明示的に示すために、ここでは、rio::をつけている。::は、パッケージ内の関数を指定するために使われる。同じ名称の関数が複数のパッケージにある場合、どのパッケージの関数を使うかを指定し、意図しないエラーを防ぐことができる。

# パッケージ読み込み
library(rio)
# データ読み込み
gravityIMFeu <- rio::import("gravityIMFeu.dta")
# データ確認
head(gravityIMFeu)

##   year iso3_o iso3_d  dist contig comlang_off   gdp_o    gdp_d gdpcap_o
## 1 1948    AUS    CAN 16025      0           1      NA       NA       NA
## 2 1949    AUS    CAN 16025      0           1 6608000 16295000    0.834
## 3 1950    AUS    CAN 16025      0           1 7042000 16979000    0.852
## 4 1951    AUS    CAN 16025      0           1 7458000 20551000    0.876
## 5 1952    AUS    CAN 16025      0           1 8007000 25117000    0.921
## 6 1953    AUS    CAN 16025      0           1 9137000 26276999    1.032
##   gdpcap_d gatt_o gatt_d wto_o wto_d eu_o eu_d rta rta_coverage rta_type
## 1       NA      1      1     0     0    0    0   0           NA       NA
## 2    1.210      1      1     0     0    0    0   0           NA       NA
## 3    1.212      1      1     0     0    0    0   0           NA       NA
## 4    1.434      1      1     0     0    0    0   0           NA       NA
## 5    1.699      1      1     0     0    0    0   0           NA       NA
## 6    1.731      1      1     0     0    0    0   0           NA       NA
##   tradeflow_imf_o tradeflow_imf_d pair treated_after treated after
## 1           26000           30100 1604             0       0     0
## 2           26600           29200 1604             0       0     0
## 3           24100           33600 1604             0       0     0
## 4           36300           48100 1604             0       0     0
## 5           19100           21100 1604             0       0     0
## 6           24600           26200 1604             0       0     0

なお、RでStata形式のデータ（.dta）を読み込むことに特化したパッケージとしては、havenがある。havenパッケージを使っても、Stata形式のデータを読み込むことができる。以下でインストールする。

install.packages("haven")

そして、havenパッケージのread_dta関数でdtaファイルを読み込む。

library(haven)
gravityIMFeu <- haven::read_dta("gravityIMFeu.dta")
head(gravityIMFeu)

3 分析の枠組み

処置群は、輸出国・輸入国の両方が、ドイツ（DEU）、フランス（FRA）、イタリア（ ITA）、オランダ（NLD）である２国間域内貿易である。制御群は、輸出国・輸入国のいずれか、もしくは両方が、ドイツ・フランス・イタリア・オランダではない２国間貿易である。1958年1月1日に、ローマ条約が発効し、欧州経済共同体（EEC、のちの欧州共同体 (EC)）発足した。そのため、1948–1957が処置前期間（before）、1958–1972が処置後期間（after）である。

	Before 1948–1957	After 1958–1972	差
EEC非加盟国	A0	A1	A1-A0
EEC加盟国	B0	B1	B1-B0
差	B0-A0	B1-A1	(B1-B0)-(A1-A0)

EEC非加盟国

A1-A0

EEC加盟国

B1-B0

差

B0-A0

B1-A1

(B1-B0)-(A1-A0)

tableにより観測数の確認を行う。

table(gravityIMFeu$treated, gravityIMFeu$after)

##    
##         0     1
##   0  7260 10890
##   1   120   180

処置群は、ドイツ（DEU）、フランス（FRA）、イタリア（ ITA）、オランダ（NLD）の4カ国の12ペアである。処置前の期間については、処置群は12x10年=120観測値、処置後の期間については、処置群は12x18年=180観測値となっている。

データを処理するために汎用的に使われているdplyrパッケージをインストールして、記述統計の確認を行う。なお、dplyrパッケージは、tidyverseパッケージの一部であるため、tidyverseパッケージをインストールしていれば、dplyrパッケージのインストールは原則不要である。

install.packages("dplyr")

まず記述統計のデータdesを作成する。
group_by関数は、データをグループ化する関数である。ここでは、treatedとafterの組み合わせでグループ化している。
以下のコードでは、gravityIMFeuデータから貿易額の処置群・制御群・処置前後のグループ別の平均値を算出するために使っている。
summarise関数は、データの要約統計量を計算する関数である。ここでは、mean関数とsd関数を使って、平均値と標準偏差を計算している。

library(dplyr)
des <- gravityIMFeu %>%
group_by(treated, after) %>%
summarise(mean = mean(log(tradeflow_imf_d)), 
          sd = sd(log(tradeflow_imf_d))) 

des # データ（処置群・制御群の貿易額の平均値と標準偏差）の表示

## # A tibble: 4 × 4
## # Groups:   treated [2]
##   treated after  mean    sd
##     <dbl> <dbl> <dbl> <dbl>
## 1       0     0  9.26 1.88 
## 2       0     1 10.1  1.84 
## 3       1     0 11.8  0.947
## 4       1     1 13.7  1.03

ここでは、tidyverseパッケージでよく用いられるパイプ演算子%>%を使っている。%>%は、左辺の結果を右辺の関数の第一引数に渡す。パイプ演算子%>%のショートカット・キーは、Macでは「Control + Shift + m」である。 Windowsでは「Shift + Alt + M 」である。なお、Base Rでは、パイプ演算子の記号は、|>である。
- 参考）「付録 B — パイプ演算子 (R)」土井翔平『R/Pythonで計量政治学入門』
パイプ演算子を用いないで書くと以下のようになる。パイプ演算子を使って作成したdesとパイプ演算子を使わないで作成したdes2は同じデータになる。

library(dplyr)
des2 <- gravityIMFeu
des2 <- group_by(gravityIMFeu, treated, after) # 第1引数は、「gravityIMFeu」
des2 <- summarise(des2, mean = mean(log(tradeflow_imf_d)), sd = sd(log(tradeflow_imf_d))) # 第1引数は、「des2」
des2 # データ（処置群・制御群の貿易額の平均値と標準偏差）の表示

## # A tibble: 4 × 4
## # Groups:   treated [2]
##   treated after  mean    sd
##     <dbl> <dbl> <dbl> <dbl>
## 1       0     0  9.26 1.88 
## 2       0     1 10.1  1.84 
## 3       1     0 11.8  0.947
## 4       1     1 13.7  1.03

処置群と制御群の貿易額の平均値をそれぞれデータtreatedとcontrolに収録する。
subset()関数は、データのサブセットを作成する関数である。ここでは、treatedが1のデータを処置群、treatedが0のデータを制御群としている。

treated <- subset(des,treated==1) # 処置群のデータ
control <- subset(des,treated==0) # 制御群のデータ

処置群と制御群の貿易の推移をplot関数で折れ線グラフにする。

plot()のデフォルトでは[0, 1]の間の0.1, 0.2,…などの数値ラベルが表示される。それを防ぐため、xaxt = "n"によりx軸の数値ラベルは非表示にしている。その上で、axis(1, 0:1)により、x軸の数値ラベルを[0, 1]のみ表示している。

# 処置群（青）の推移。
plot(treated$after, treated$mean, 
     type = "o", ylim=c(8,15), col=4, 
     xlab = "period", ylab="Trade values",  xaxt = "n") 

# x軸の数値ラベルを[0, 1]にしている。
axis(1, 0:1) # axisの最初の引数1は、x軸を意味する。

# 制御群(黒)の推移。
points(control$after,control$mean,type = "o", col=1) 


text(0.5,13,"treated") # 処置群のラベル
text(0.5,10,"control") # 制御群のラベル

text(0, 11.8+0.5, "11.8") # 処置群のbeforeの平均値
text(1, 13.7+0.5, "13.7") # 処置群のafterの平均値

text(0, 9.3+0.5, "9.3") # 制御群のbeforeの平均値
text(1, 10.1+0.5, "10.1") # 制御群のafterの平均値

4 伝統的重力方程式に基づく推定

4.1 単純な2 x 2 DiD – `Base R`

まず、OLS推定で、単純な2 x 2 の古典的なDiDで、EU加盟の効果を推定する。処置変数treatedと処置後を表すafterの交差項treated * afterは、処置後の期間のEEC域内貿易であれば1をとる。

did <- lm(log(tradeflow_imf_d) ~  treated * after
          + log(gdp_o) + log(gdp_d) + log(dist), 
          data = gravityIMFeu)

4.2 TWFE (OLS) – `Base R`

観測単位（貿易国ペア）と時間（年）両方の固定効果（ダミー変数）を加えて、処置効果を推定する「２元（双方向）固定効果モデル」 (TWFE model: two-way fixed effects model) が近年使用されることが多い。

TWFEモデルをBase RでOLS推定するコードは以下の通りである。推定にはかなり時間がかかる。

twfe <- lm(log(tradeflow_imf_d) ~ treated_after 
           + log(gdp_o) + log(gdp_d)
           + factor(pair) + factor(year), data = gravityIMFeu)

4.3 TWFE (OLS) – `fixest`

TWFEモデルをfixestでOLS推定するコードは以下の通りである。なお、fixestパッケージを用いていることを明示するために、fixest::を前置する。本来は、library(fixest)でパッケージを読み込めば、fixest::を前置する必要はない。

twfe1 <- fixest::feols(log(tradeflow_imf_d) ~ treated_after 
           + log(gdp_o) + log(gdp_d)
           | pair + year, 
           vcov = ~ iso3_o + iso3_d,
           data = gravityIMFeu)

4.4 TWFE (PPML) –`fixest`

TWFEモデルをfixestでPPML推定するコードは以下の通りである。

twfe2 <- fixest::fepois(tradeflow_imf_d  ~ treated_after 
                + log(gdp_o) + log(gdp_d) 
                | pair + year, 
                vcov = ~ iso3_o + iso3_d,
                data = gravityIMFeu)

5 固定効果アプローチに基づく推定

5.1 TWFE (PPML) –`fixest`

TWFEモデルをfixestで固定効果アプローチでPPML推定するコードは以下の通りである。

twfe3 <- fixest::fepois(tradeflow_imf_d  ~ treated_after 
                | iso3_o^year + iso3_d^year + pair + year,
                vcov = ~ iso3_o + iso3_d,
                data = gravityIMFeu)

ここで、以下のようにyearを固定効果に指定しなくても同じ結果が得られる。これは、iso3_o^yearもしくは、iso3_d^yearを固定効果に指定しているため、yearの固定効果も自動的に含まれるからである。ただ、推定結果表で、yearの固定効果を明示するために、上のコードでは、yearを固定効果に指定した。

twfe3 <- fixest::fepois(tradeflow_imf_d  ~ treated_after 
                | iso3_o^year + iso3_d^year + pair,
                vcov = ~ iso3_o + iso3_d,
                data = gravityIMFeu)

6 結果の比較

modelsummaryを使って、推定結果の比較を行う。
coef_renameで係数の表示名を統一化している。
そのうえで、kableExtraパッケージを使って、表の見出しや脚注を追加し、表の見栄えを整えている。
パッケージのインストールは以下のコードで行える。バージョンによって、modelsummaryやkableExtraの機能が異なるため、バージョンを指定してインストールする。

# 古いパッケージのインストールのため
install.packages("remotes")

# modelsummaryパッケージをインストールする。
remotes::install_version(package = "modelsummary", version = "1.4.3", 
                         dependencies = FALSE,
                         repos = "http://cran.us.r-project.org")

# LaTexへの表の出力のために、kableExtraパッケージをインストールする。
#install.packages("kableExtra")
remotes::install_version(package = "kableExtra", version = "1.4.0", 
                         dependencies = FALSE,
                         repos = "http://cran.us.r-project.org")

library(modelsummary)
library(kableExtra)

# longnote
longnote <- "Gravity equation estimates are reported. Standard errors are clustered at country-level and presented in parentheses."


coef_rename = c("treated × after" = "treated × after", 
                "treated:after" = "treated × after", 
                "treated_after" = "treated × after")

models <- list("2x2 DiD \ (GDP)" = did, 
               "TWFE-OLS \ (GDP)" = twfe1,
               "TWFE-PPML \ (GDP)" = twfe2, 
               "TWFE-PPML \ (Fixed Effects)" = twfe3)

modelsummary::modelsummary(models,
             stars = TRUE,
             coef_omit = "Intercept|.*factor",
             coef_rename = coef_rename,
             gof_omit = 'DF|Deviance|Log.Lik.|RMSE|AIC|BIC|R2 Within|R2 Within Adj. ') %>%
kableExtra::add_header_above(c(" " = 1, 
                               "OLS" = 2,  
                               "PPML" = 2)) %>%
kableExtra::add_header_above(c(" " = 1, 
                               "(1)" = 1,  
                               "(2)" = 1, 
                               "(3)" = 1, 
                               "(4)" = 1)) %>% 
kableExtra::footnote(general = longnote, 
                     threeparttable = TRUE)

	(1)	(2)	(3)	(4)
	OLS		PPML
	2x2 DiD (GDP)	TWFE-OLS (GDP)	TWFE-PPML (GDP)	TWFE-PPML (Fixed Effects)
treated	0.001
	(0.129)
after	−0.533***
	(0.019)
log(gdp_o)	0.723***	0.564***	0.859***
	(0.005)	(0.156)	(0.105)
log(gdp_d)	0.720***	0.532***	0.622***
	(0.005)	(0.138)	(0.097)
log(dist)	−0.505***
	(0.009)
treated × after	0.919***	0.792***	0.699***	0.504***
	(0.151)	(0.068)	(0.112)	(0.144)
Num.Obs.	16319	16319	16319	18450
R2	0.688	0.935	0.979
R2 Adj.	0.688	0.932	0.979
F	5999.137
Std.Errors		by: iso3_o & iso3_d	by: iso3_o & iso3_d	by: iso3_o & iso3_d
FE: pair		X	X	X
FE: year		X	X	X
FE: iso3_o^year				X
FE: iso3_d^year				X
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001
Note:
Gravity equation estimates are reported. Standard errors are clustered at country-level and presented in parentheses.

modelplotを用いて係数プロットを作成し、モデルを比較する。ここで、coef_renameを使って、変数の表示名を統一化している。また、coef_omit = "^(?!.*treated|.*after)"によって、treatedもしくはafterを含まない変数はのぞいている。

coef_rename = c("treated × after" = "treated × after", 
                "treated:after" = "treated × after", 
                "treated_after" = "treated × after")

modelsummary::modelplot(models, 
          coef_omit = "^(?!.*treated|.*after)" ,
          coef_rename =coef_rename)

参考文献

Anderson, James E, and Eric Van Wincoop. 2003. “Gravity with Gravitas: A Solution to the Border Puzzle.” American Economic Review 93 (1): 170–92. https://doi.org/10.1257/000282803321455214.

Berge, Laurent, and Grant McDermott. 2023. “Fast Fixed-Effects Estimation: Short Introduction.” https://cran.r-project.org/web/packages/fixest/vignettes/fixest_walkthrough.html.

Goodman-Bacon, Andrew. 2021. “Difference-in-Differences with Variation in Treatment Timing.” Journal of Econometrics 225 (2): 254–77.

Roth, Jonathan, Pedro HC Sant’Anna, Alyssa Bilinski, and John Poe. 2023. “What’s Trending in Difference-in-Differences? A Synthesis of the Recent Econometrics Literature.” Journal of Econometrics 235 (2): 2218–44.

補論1

推定結果のLaTexへの出力を行う。

推定結果表は、基本的にmodelsummaryパッケージで作成。
output = "latex"でLaTexへの出力を指定。
save_kable("table1.tex")で、LaTexファイルを保存する。

modelsummary::modelsummary(models, 
             output = "latex",
             stars = TRUE,
             coef_omit = "Intercept|.*factor",
             coef_rename = coef_rename,
             gof_omit = 'DF|Deviance|Log.Lik.|RMSE|AIC|BIC|R2 Within|R2 Within Adj. ') %>%
kableExtra::add_header_above(c(" " = 1, 
                               "OLS" = 2,  
                               "PPML" = 2)) %>%
kableExtra::add_header_above(c(" " = 1, 
                               "(1)" = 1,  
                               "(2)" = 1, 
                               "(3)" = 1, 
                               "(4)" = 1)) %>% 
kableExtra::footnote(general = longnote, 
                     threeparttable = TRUE) %>% 
kableExtra::save_kable("table1.tex")

補論2

モデルの名前を推定結果表の表頭に記す別の方法として、tibbleパッケージのtibble関数で、表頭のデータフレームを作成することもできる。この場合、attr(rows, 'position') <- c(1, 1)と、一番最初にモデルの名前を表示するように、設定する。

tibbleをインストールしていない場合は、以下でインストールする。

install.packages("tibble")

以下のコードで、モデルの名前を推定結果表の表頭に記すことができる。

library(modelsummary)
library(tibble)
rows <- tibble(
  a = "Model",
  b = "2x2 DiD",
  c = "TWFE-OLS",
  d = "TWFE-PPML",
  e = "TWFE-PPML",
)
attr(rows, 'position') <- c(1, 1)
modelsummary(list(did, twfe1, twfe2, twfe3), 
             stars = TRUE,
             coef_omit = "Intercept|.*factor",
             gof_omit = 'DF|Deviance|Log.Lik.|RMSE|AIC|BIC|R2 Within|R2 Within Adj. ',
             add_rows = rows)

	(1)	(2)	(3)	(4)
Model	2x2 DiD	TWFE-OLS	TWFE-PPML	TWFE-PPML
treated	0.001
	(0.129)
after	−0.533***
	(0.019)
log(gdp_o)	0.723***	0.564***	0.859***
	(0.005)	(0.156)	(0.105)
log(gdp_d)	0.720***	0.532***	0.622***
	(0.005)	(0.138)	(0.097)
log(dist)	−0.505***
	(0.009)
treated × after	0.919***
	(0.151)
treated_after		0.792***	0.699***	0.504***
		(0.068)	(0.112)	(0.144)
Num.Obs.	16319	16319	16319	18450
R2	0.688	0.935	0.979
R2 Adj.	0.688	0.932	0.979
F	5999.137
Std.Errors		by: iso3_o & iso3_d	by: iso3_o & iso3_d	by: iso3_o & iso3_d
FE: pair		X	X	X
FE: year		X	X	X
FE: iso3_o^year				X
FE: iso3_d^year				X
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001

重力方程式を用いた2元固定効果差の差推定

田中鮎夢

2024-03-14

1 はじめに

2 貿易データ

3 分析の枠組み

4 伝統的重力方程式に基づく推定

4.1 単純な2 x 2 DiD – `Base R`

4.2 TWFE (OLS) – `Base R`

4.3 TWFE (OLS) – `fixest`

4.4 TWFE (PPML) –`fixest`

5 固定効果アプローチに基づく推定

5.1 TWFE (PPML) –`fixest`

6 結果の比較

参考文献

補論1

補論2

重力方程式を用いた2元固定効果差の差推定

田中 鮎夢

2024-03-14

1 はじめに

2 貿易データ

3 分析の枠組み

4 伝統的重力方程式に基づく推定

4.1 単純な2 x 2 DiD – Base R

4.2 TWFE (OLS) – Base R

4.3 TWFE (OLS) – fixest

4.4 TWFE (PPML) –fixest

5 固定効果アプローチに基づく推定

5.1 TWFE (PPML) –fixest

6 結果の比較

参考文献

補論1

補論2

田中鮎夢

4.1 単純な2 x 2 DiD – `Base R`

4.2 TWFE (OLS) – `Base R`

4.3 TWFE (OLS) – `fixest`

4.4 TWFE (PPML) –`fixest`

5.1 TWFE (PPML) –`fixest`