fixestThis version: 2024-06-04
Anderson and Van Wincoop (2003) が指摘した多角的貿易抵抗指数（物価効果）を制御するために、1時点のクロスセクションの貿易データであれば、輸出国のGDPの代わりに輸出国固定効果、輸入国のGDPの代わりに輸入国固定効果を説明変数に加えれば良い。２時点以上のパネルの貿易データであれば、毎年GDPや物価は変化するので、GDPの代わりに輸出国（輸入国）固定効果と年次固定効果の交差項を説明変数に加えることで、時変の輸出国（輸入国）属性を制御する。このように、GDPの代わりに、輸出国（輸入国）固定効果を用いるアプローチを「固定効果アプローチ」と呼んでいる。
さらに、重力方程式に、貿易ペアの固定効果を含めることがある。例えば、イギリスとアメリカという貿易国ペア、日本とアメリカという貿易国ペアごとにダミー変数を作成して、重力方程式の説明変数に加えるということである。これにより、貿易ペア固有の時間不変の要因を制御できる。貿易ペアの固定効果を推定に加えることは、計量経済学のテキストで「固定効果法（固定効果モデル）」と呼ばれているパネルデータ推定法と同じである。ここでいう「固定効果法」は、先に述べた多角的貿易抵抗指数を制御するための「固定効果アプローチ」とは意味が違うことに注意が必要である。
このように、時変の輸出国（輸入国）固定効果、貿易ペアの固定効果、さらには年次固定効果を加えて、重力方程式をポワソン擬似最尤（PPML）で推定するのが、現代の重力方程式の標準的な推定方法である。多数の固定効果(ダミー変数)を含むため、Base
Rでは計算上困難であるため、本ページでは、Berge
and McDermott (2023)
によって開発されたfixestパッケージを用いる。
パッケージのインストールは以下のコードで行われる。
install.packages("fixest")パッケージの読み込みを行う。
# パッケージの読み込み
library(fixest)gravity_rta.csvには、フランスの研究機関CEPIIが作成している、Gravityという無料のデータベースから、２国間の貿易額やRTAの有無など必要最小限のデータを抜き出している。gravity_rta.csvには、2010–2018年の9年間の世界の2国間の貿易額が収録されている。Gravityデータベースは定期的に更新されている。本ページで用いるデータは、「Gravity_dta_V202211」版であり、2024年3月に入手した。
Rは大規模なデータの読み込みには非常に時間がかかる。事実上読み込めないことがある。今回使うデータは、readxlパッケージのread_excelではexcel形式で読み込めない。そのため、Dowle et al. (2019)
が開発した大規模なデータの読み込みに対応しているdata.tableパッケージをまずインストールする。
install.packages("data.table")そして、data.tableパッケージのfread関数でcsvファイルを読み込む。
library(data.table)
gravity_rta <- fread("gravity_rta.csv")
head(gravity_rta)## year iso3_o iso3_d rta tradeflow_baci pair
## <int> <char> <char> <int> <num> <int>
## 1: 2017 ABW AFG 0 160.164 2
## 2: 2018 ABW AFG 0 1023.860 2
## 3: 2010 ABW AGO 0 50257.141 3
## 4: 2011 ABW AGO 0 0.142 3
## 5: 2012 ABW AGO 0 0.329 3
## 6: 2013 ABW AGO 0 25.000 3我々は、貿易に対するRTA（regional trade agreement）の効果を検証するために、重力モデルを推定することにする。従属変数は二国間の貿易水準であり、独立変数は二国間のRTAダミーである。以下のように定式化する：
\[ E\left(\text { Trade }_{i, j, t}\right)=\gamma_{i,t}^{\text {Exporter }} \times \gamma_{j,t}^{\text {Importer }} \times \gamma_{t}^{\text {Year }} \times \gamma_{i,j}^{\text {Pair}} \times \text { RTA }_{i j, t}^{\beta} \] ここで、添え字\(i\)、\(j\)、\(t\)はそれぞれ輸出国、輸入国、年を表し、\(\gamma\)はこれらのグループの固定効果である。ここで\(\beta\)は注目される処置効果である。
Silva and Tenreyro (2006) に従い、ポワソン擬似最尤（PPML）推定を用いて重力方程式を推定することを考える。これは次の関係を導く: \[ E\left(\text { Trade }_{i, j, t}\right)=\exp \left(\gamma_{i,t}^{\text {Exporter }}+\gamma_{j,t}^{\text {Importer }}+\gamma_{t}^{\text {Year }}+ \gamma_{i,j}^{\text {Pair}} + \beta \times \text { RTA }_{i j, t}\right) \]
まずは、OLSで推定する。線形で推定するために、変数を対数にする必要がある：
gravity_ols = feols(log(tradeflow_baci) ~ rta | iso3_o^year + iso3_d^year + pair, vcov = ~ iso3_o + iso3_d, data = gravity_rta)結果はprintもしくはsummaryで直接表示できる：
print(gravity_ols)## OLS estimation, Dep. Var.: log(tradeflow_baci)
## Observations: 282,673
## Fixed-effects: iso3_o^year: 2,017, iso3_d^year: 2,017, pair: 40,237
## Standard-errors: Clustered (iso3_o & iso3_d)
## Estimate Std. Error t value Pr(>|t|)
## rta 0.040042 0.044688 0.896029 0.37119
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 1.32381 Adj. R2: 0.894409
## Within R2: 5.43e-6Stataで同じ推定を行うにはreghdfeをインストールする。
ssc install reghdfeデータを読み込んだ上で、以下のコードを実行して、対数をとる。
import delimited using https://ayumu-tanaka.github.io/teaching/gravity_rta.csv,clear
g lntradeflow_baci=ln(tradeflow_baci)さらに、文字変数をabsorbで指定することができないので、輸出国（輸入国）ISOコードを数値変数に変換しておく。
encode iso3_o,gen(origin)
encode iso3_d,gen(destination)最小二乗法の実行は、以下の通りである。
reghdfe lntradeflow_baci rta, absorb(i.origin##i.year i.destination##i.year i.pair) vce(cluster origin destination)Stataでは、A##Bは、Aの固定効果、Bの固定効果、AとBの交差項の固定効果全てを使うことを意味する。
ポアソン尤度を用いたモデルの推定は以下の通りである：
gravity_pois = fepois(tradeflow_baci ~ rta | iso3_o^year + iso3_d^year + pair, vcov = ~ iso3_o + iso3_d, data = gravity_rta)結果はprintもしくはsummaryで直接表示できる：
print(gravity_pois)## Poisson estimation, Dep. Var.: tradeflow_baci
## Observations: 282,673
## Fixed-effects: iso3_o^year: 2,017, iso3_d^year: 2,017, pair: 40,237
## Standard-errors: Clustered (iso3_o & iso3_d)
## Estimate Std. Error z value Pr(>|z|)
## rta 0.078878 0.020277 3.89005 0.00010022 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Log-Likelihood: -2.797e+9 Adj. Pseudo R2: 0.99451
## BIC: 5.594e+9 Squared Cor.: 0.995649printは，係数推定値と標準誤差，およびいくつかの他の情報を報告する。適合度の情報のうち、二乗相関は従属変数と予測変数の相関に対応し、OLS推定におけるR2乗の考え方を反映している。
Stataで同じ推定を行うには次のコードをppmlhdfeをインストールする。
ssc install ppmlhdfeデータを読み込み、以下のコードで、PPML推定を行う。
import delimited using https://ayumu-tanaka.github.io/teaching/gravity_rta.csv,clear
encode iso3_o,gen(origin)
encode iso3_d,gen(destination)
ppmlhdfe tradeflow_baci rta,absorb(i.origin##i.year i.destination##i.year i.year i.pair) vce(cluster origin destination)ここで、いくつかの推定結果をコンパクトに概観するため、関数etable
を使う。この関数は、複数の固定効果推定の結果をdata.frameに要約する。推定結果をまとめて見るには、次のようにタイプするだけでよい：
etable(gravity_ols, gravity_pois,
vcov = ~ iso3_o + iso3_d,
headers = c("OLS","PPML"))## gravity_ols gravity_pois
## OLS PPML
## Dependent Var.: log(tradeflow_baci) tradeflow_baci
##
## rta 0.0400 (0.0447) 0.0789*** (0.0203)
## Fixed-Effects: ------------------- ------------------
## iso3_o-year Yes Yes
## iso3_d-year Yes Yes
## pair Yes Yes
## _______________ ___________________ __________________
## Family OLS Poisson
## S.E.: Clustered by: iso3_o & iso3_d by: iso3_o & iso3_d
## Observations 282,673 282,673
## Squared Cor. 0.91095 0.99565
## Pseudo R2 0.41574 0.99451
## BIC 1,516,456.0 -2,147,483,648.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1PPML推定の結果得られるRTAの推定係数は、0.0789であり、統計的に有意である。RTAは0か1のダミー変数であるので、RTAが貿易を1.082倍(\(=\exp(\beta \times RTA)=\exp(0.0789 \times 1)\)倍)程度増加させる効果があることを示している。Rでは、以下のように、RTAの効果を計算できる。
exp(0.0789)## [1] 1.082096これまで、複数の推定結果をRコンソールで報告する方法を見てきた。ここで、同じ関数
etable を使用して、結果を LaTex
形式の表にエクスポートすることができる。
etable(gravity_ols, gravity_pois, cluster = ~iso3_o + iso3_d, file = "04_Tables2.tex", replace = TRUE)この例では、2つの推定結果を含む1つの表がLaTexの表に直接エクスポートされ、“Tables2.tex”というファイルに格納されている。引数ファイルが存在する場合、LaTexフォーマットがデフォルトになるため、引数tex=TRUEを使用する必要がない。このファイルは、replace=TRUEという引数のおかげで、毎回再作成される。