Gravitationsgleichungen schätzen

11.05.2020

Wie schätzt man Gravitationsgleichungen richtig? Und wie interpretiert man Gravity-Koeffizienten? In dieser Einheit dreht sich alles um unterschiedliche Schätzmethoden.

Vorlesung

Fullscreen PDF Audio

Anwendung

Wir nutzen nun das in der Vorlesung erarbeitete verallgemeinerte strukturelle Gravitationsmodell zur Simulation der Effekte von Handelspolitik, am Beispiel des Brexits.

Daten laden und vorbereiten

Starten Sie RStudio und legen Sie neues Projekt an. Wie in der letzten Einheit benötigen wir eine Reihe von Packages.

pacman::p_load(readr)
pacman::p_load(data.table)
pacman::p_load(stringr)
pacman::p_load(ggplot2)
pacman::p_load(scales)
pacman::p_load(lfe)
pacman::p_load(alpaca)

Diesmal werden die Packages zwar wieder mit p_load aus dem Paket pacman geladen, allerdings ohne dieses erst komplett zu laden. Wir nutzen lediglich die p_load Funktion mittels dem Operator ::.

Wir nutzen heute einen anderen Datensatz als zuvor, nämlich einen, der auch eine Zeitdimension hat:

trade_gravity.rds: Bilaterale Handelsdaten und Gravity-Variablen vom CEPII, ursprünglich aus Head, Mayer & Ries (2011)¹

Laden Sie die Daten mittels read_rds():

# Daten laden ####
data_gravity = read_rds("input/trade_gravity.rds")

Nebenbemerkung: Zur Übersicht setze ich in meinem R-Code häufig section marks mittels # Titel ####.²

Neue Daten kennenlernen

Back to work. Machen Sie sich ein wenig mit den Daten vertraut, in dem sie mit head(data_gravity) sich die Variablen anschauen, und vielleicht auch einige Zusammenhänge mit ggplot() plotten. Z.b. die Anzahl der Nullen über die Zeit, wie in der Vorlesung gesehen:

# Nullen über die Zeit (Extensiver Rand)
plot_data = data_gravity[, .(share_zero = sum(value == 0) / .N), by = year]
plot = ggplot(plot_data) +
  theme_minimal() +
  geom_line(aes(x = year, y = share_zero)) +
  scale_x_continuous() +
  scale_y_continuous("Anteil Länderpaare ohne Handel", labels = scales::percent_format(accuracy = 1)) +
  theme(axis.title.x = element_blank())
print(plot)
ggsave(plot, filename = "output/nullen.png", width = 20, height = 10, units = "cm")
rm(plot, plot_data)

… oder das Gesamtvolumen des Welthandels über die Zeit:

# Intensiver Rand
plot_data = data_gravity[, .(value = sum(value)), by = year]
plot = ggplot(plot_data) +
  theme_minimal() +
  geom_line(aes(x = year, y = value / 1000)) +
  scale_x_continuous() +
  scale_y_continuous("Weltweites Handelsvolumen in Mrd. USD") +
  theme(axis.title.x = element_blank())
print(plot)
ggsave(plot, filename = "output/handel_linear.png", width = 20, height = 10, units = "cm")

bzw. besser in mit logarithmierter y-Achse:

plot = ggplot(plot_data) +
  theme_minimal() +
  geom_line(aes(x = year, y = value / 1000)) +
  scale_x_continuous() +
  scale_y_log10("Weltweites Handelsvolumen in Mrd. USD") +
  theme(axis.title.x = element_blank())
print(plot)
ggsave(plot, filename = "output/handel_log10.png", width = 20, height = 10, units = "cm")
rm(plot, plot_data)

OLS ohne multilaterale Resistanzterme

Nun aber zur ersten Regression. Wir schätzen jetzt, wie in der zweiten Vorlesung, ein naives Gravitationsmodell. Da es, wie wir oben sehen, eine Menge Beobachtungen mit value == 0 gibt, und wir aber mit log(value) schätzen, schließen wir diese aus der Regression aus. Das sieht dann so aus:

# 1. OLS ohne multilaterale Resistanzterme ####
reg1 = lm(log(value) ~
            log(gdp_o) + log(gdp_d) +
            log(distw) + contig +
            colony + comlang_off +
            fta + comcur,
          data = data_gravity[value > 0])
summary(reg1)

Der Regressionsoutput sollte in etwa so aussehen:

Call:
lm(formula = log(value) ~ log(gdp_o) + log(gdp_d) + log(distw) +
    contig + colony + comlang_off + fta + comcur, data = data_gravity[value >
    0])

Residuals:
    Min      1Q  Median      3Q     Max
-21.140  -1.317   0.296   1.617   9.590

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -5.792602   0.040090 -144.49   <2e-16 ***
log(gdp_o)   0.871607   0.001404  620.88   <2e-16 ***
log(gdp_d)   0.704754   0.001338  526.89   <2e-16 ***
log(distw)  -1.029101   0.004297 -239.47   <2e-16 ***
contig       0.681060   0.019102   35.65   <2e-16 ***
colony       1.786580   0.019920   89.69   <2e-16 ***
comlang_off  0.392085   0.008532   45.95   <2e-16 ***
fta          0.563299   0.015527   36.28   <2e-16 ***
comcur       0.759907   0.026123   29.09   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.405 on 624136 degrees of freedom
  (85428 observations deleted due to missingness)
Multiple R-squared:  0.5229,	Adjusted R-squared:  0.5229
F-statistic: 8.55e+04 on 8 and 624136 DF,  p-value: < 2.2e-16

Die Koeffizienten sind alle ungefähr so wie erwartet. Ein wenig erstaunlich ist vielleicht, dass die Koeffizienten für log(gdp_o) und log(gdp_d) nicht so nah an 1 sind.³

OLS mit fixen Exporteur- und Importeur-spezifischen Effekten

Multilaterale Resistanzterme sind aber wichtig. In der Vorlesung haben wir gesehen, dass wenn man sich allein für den Einfluss von bilateralen Variablen interessiert — was häufig der Fall ist — dann sind fixe Effekte das Mittel der Wahl. Da wir eine Zeitdimension haben, benötigen wir fixe Effekte für Land-Jahr Kombinationen, die man ganz einfach so generieren kann:

# 2. OLS mit fixen Exporteur- und Importeur-spezifischen Effekten ####
data_gravity[, iso_o_year := str_c(iso_o, year)]
data_gravity[, iso_d_year := str_c(iso_d, year)]

In den Daten haben wir nun also zwei weitere Variablen, die einfach der zusammengesetzte String aus der iso_o, bzw. iso_d, und year Variablen sind:

> head(data_gravity)
   year iso_o iso_d ... iso_o_year iso_d_year
1: 1999   ABW   AGO ...    ABW1999    AGO1999
2: 2000   ABW   AGO ...    ABW2000    AGO2000
3: 2001   ABW   AGO ...    ABW2001    AGO2001
4: 2002   ABW   AGO ...    ABW2002    AGO2002
5: 2003   ABW   AGO ...    ABW2003    AGO2003
6: 2004   ABW   AGO ...    ABW2004    AGO2004

Zur schnellen Schätzung eines solchen “fixed effects”-Modells am besten die felm() Funktion aus dem lfe Package nehmen. Die Syntax ist sehr ähnlich der von lm, nur muss man die fixen Effekte in die Schätzgleichung mittels “|”-Operator einfügen. Das sieht dann so aus:

reg2 = felm(log(value) ~
              log(distw) + contig + colony + comlang_off + fta + comcur |
              iso_o_year + iso_d_year,
            data = data_gravity[value > 0])
summary(reg2)

Das Ergebnis sollte so aussehen:

Call:
   felm(formula = log(value) ~ log(distw) + contig + colony + comlang_off +
         fta + comcur | iso_o_year + iso_d_year, data = data_gravity[value >      0])

Residuals:
     Min       1Q   Median       3Q      Max
-20.7024  -1.0019   0.1105   1.1174  12.7477

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
log(distw)  -1.269274   0.003812 -332.98   <2e-16 ***
contig       0.444242   0.015036   29.54   <2e-16 ***
colony       1.265863   0.016176   78.26   <2e-16 ***
comlang_off  0.560282   0.007543   74.28   <2e-16 ***
fta          0.606633   0.013277   45.69   <2e-16 ***
comcur       0.941889   0.019482   48.35   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.898 on 690130 degrees of freedom
Multiple R-squared(full model): 0.7073   Adjusted R-squared: 0.699
Multiple R-squared(proj model): 0.2558   Adjusted R-squared: 0.2348
F-statistic(full model):85.77 on 19442 and 690130 DF, p-value: < 2.2e-16
F-statistic(proj model): 3.954e+04 on 6 and 690130 DF, p-value: < 2.2e-16

Die geschätzten Koeffizienten haben sich zum Teil deutlich verändert. Das sollte kein Wunder sein, denn die fixen Effekte absorbieren ja nun alles Exporteur- und Importeur-spezifisches: Sowohl wirtschaftliche Größe, als auch die wichtigen multilateralen Resistanzterme, und alle möglichen anderen Länder-spezifischen Charakteristika.

OLS mit fixen Exporteur-, Importeur-spezifischen, und bilateralen Effekten

In der Vorlesung hatten wir gesehen, dass zusätzliche bilaterale fixe Effekte helfen, unbeobachtbare bilaterale Charakteristika zu absorbieren. Wir generieren also eine weitere zusätzliche Variable iso_o_iso_d und ändern die Schätzgleichung entsprechend:

# 3. OLS mit fixen Exporteur-, Importeur-spezifischen, und bilateralen Effekten ####
data_gravity[, iso_o_iso_d := str_c(iso_o, iso_d)]
reg3 = felm(log(value) ~
              fta + comcur |
              iso_o_year + iso_d_year + iso_o_iso_d,
            data = data_gravity[value > 0])
summary(reg3)

Das Ergebnis sieht so aus:

Call:
   felm(formula = log(value) ~ fta + comcur | iso_o_year + iso_d_year +
         iso_o_iso_d, data = data_gravity[value > 0])

Residuals:
     Min       1Q   Median       3Q      Max
-18.0374  -0.5939   0.0445   0.6849  10.2628

Coefficients:
       Estimate Std. Error t value Pr(>|t|)    
fta     0.45483    0.01355   33.56   <2e-16 ***
comcur  0.30739    0.02410   12.75   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.379 on 658640 degrees of freedom
Multiple R-squared(full model): 0.8525   Adjusted R-squared: 0.8411
Multiple R-squared(proj model): 0.001971   Adjusted R-squared: -0.07521
F-statistic(full model):74.77 on 50932 and 658640 DF, p-value: < 2.2e-16
F-statistic(proj model): 650.4 on 2 and 658640 DF, p-value: < 2.2e-16
*** Standard errors may be too high due to more than 2 groups and exactDOF=FALSE

Nur noch Koeffizienten für zeitvariante Variablen können jetzt geschätzt werden. Invariante Variablen, wie z.B. die bilaterale Distanz, sind kollinear zu den bilateralen fixen Effekten. Die geschätzten Koeffizienten haben sich im Vergleich zur vorherigen Schätzung deutlich geändert: Der geschätzte Koeffizient für den Effekt eines Freihandelsabkommens fta ist um ca. 25 % kleiner als vorher. Der Einfluss für eine gemeinsame Währung reduziert sich um fast 70 %! Gleichzeitig sehen wir, dass der $R^2$ deutlich angestiegen ist, auf über 0.85, von vormals 0.71, bzw. 0.52.

PPML mit fixen Exporteur-, Importeur-spezifischen, und bilateralen Effekten

Die de-facto Standartmethode zur Schätzung von Gravitationsmodellen ist ein Poisson Pseudo Maximum Likelihood Schätzer, kurz PPML. Die Schätzung lässt sich leicht implementieren mit dem alpaca Paket und der feglm Funktion. Lediglich die family Option muss man noch setzen, und zwar auf family = poisson().⁴

# 4. PPML mit fixen Exporteur-, Importeur-spezifischen, und bilateralen Effekten ####
reg4 = feglm(value ~
               fta + comcur |
               iso_o_year + iso_d_year + iso_o_iso_d,
             family = poisson(),
             data = data_gravity)
summary(reg4)

Das Ergebnis sieht so aus:

poisson - log link

value ~ fta + comcur | iso_o_year + iso_d_year + iso_o_iso_d

Estimates:
        Estimate Std. error z value Pr(> |z|)    
fta    0.2985891  0.0006655  448.64    <2e-16 ***
comcur 0.0109309  0.0008074   13.54    <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

residual deviance= 14580343.80,
null deviance= 1204872429.60,
n= 1112930, l= [9726, 9770, 31495]

( 91741 observation(s) deleted due to perfect classification )

Number of Fisher Scoring Iterations: 13

Die geschätzten Koeffizienten sind nun erheblich kleiner! Der Einfluss einer gemeinsamen Währung ist zwar statistisch signifikant unterschiedlich von Null, aber ökonomisch kaum relevant.⁵

Vergleich der Vorhersage

Am Ende der letzten Anwendung tat sich die Frage auf, wie gut die Vorhersagemöglichkeit des Modells überhaupt ist. Diese war nicht sonderlich gut. Vergleichen wir jetzt noch einmal die tatsächlich beobachteten Handelsströme mit den Vorhersagen der obigen Modelle. Dazu fügen wir die “fitted values” der Modelle als neue Variablen den Daten hinzu:

# Wie gut ist die Vorhersage? ####
data_gravity[value > 0 & !is.na(gdp_o) & !is.na(gdp_d), fitted_value_reg1 := exp(fitted(reg1))]
data_gravity[value > 0, fitted_value_reg2 := exp(fitted(reg2))]
data_gravity[value > 0, fitted_value_reg3 := exp(fitted(reg3))]
data_gravity = merge(data_gravity,
                     reg4$data[, .(iso_o_year, iso_d_year, iso_o_iso_d, fitted_value_reg4 = fitted(reg4))],
                     by = c("iso_o_year", "iso_d_year", "iso_o_iso_d"),
                     all.x = T)

Für das mit alpaca geschätzte Modell ist das leider (noch) nicht so einfach wie für die anderen beiden Modelle.

Nehmen wir nun das Beispiel USA im Jahr 2000:⁶

p_load(wesanderson)
top_ten_usa_2000 = data_gravity[year == 2000 & iso_o == "USA", .(value = sum(value)), by = iso_d][order(-value)][1:10][,iso_d]
data_plot = data_gravity[year == 2000 & iso_o == "USA" & iso_d %in% top_ten_usa_2000,
                         .(iso_d, value,
                           fitted_value_reg1,
                           fitted_value_reg2,
                           fitted_value_reg3,
                           fitted_value_reg4)]
data_plot = melt(data_plot, id.vars = "iso_d")
data_plot[, name_land := countrycode(iso_d, "iso3c", "country.name.de")]

plot = ggplot(data_plot) +
  theme_minimal() +
  geom_bar(aes(x = reorder(name_land, -value), y = value, group = variable, fill = variable),
           stat = "identity", position = "dodge") +
  scale_y_continuous("Exportvolumen", labels = dollar_format(accuracy = 1)) +
  scale_x_discrete(NULL) +
  scale_fill_manual(NULL, labels = c("Beobachtete Exporte",
                                     "Geschätzte Exporte (Naiv)",
                                     "Geschätzte Exporte (Imp. & Exp. FE)",
                                     "Geschätzte Exporte (Imp., Exp. & Bil. FE)",
                                     "Geschätzte Exporte (PPML mit FE)"),
                    values = wes_palette(n = 5, "Darjeeling1")) +
  guides(fill = guide_legend(nrow = 2)) +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 45, hjust = 1))
print(plot)
ggsave(plot, filename = "output/geschaetzte_beobachtete_exporte_usa_2000_top10.png", width = 20, height = 10, units = "cm")
rm(data_plot, plot)

Die gelben Balken (naives Gravitationsmodell) sind, wie in der letzten Woche, weiterhin nicht sonderlich ähnlich den roten Balken (tatsächlichen Handelsströme). Die hellblauen Balken (PPML Schätzung) sind aber super!

Noch einmal der direkte Vergleich von naiver Gravity-Schätzung und PPML:

data_plot = data_gravity[iso_o == "USA" & year == 2000, .(value, fitted_value_reg1, fitted_value_reg4)]
plot = ggplot(data_plot) +
  theme_minimal() +
  geom_point(aes(x = value, y = fitted_value_reg1, color = "Geschätzte Exporte (Naiv)")) +
  geom_point(aes(x = value, y = fitted_value_reg4, color = "Geschätzte Exporte (PPML mit FE)")) +
  scale_x_log10("Beobachtete Exporte") +
  scale_y_log10("Geschätzte Exporte") +
  scale_color_manual(NULL, values = wes_palette(n = 2, "GrandBudapest1"))
print(plot)
ggsave(plot, filename = "output/geschaetzte_beobachtete_exporte_usa_2000.png", width = 20, height = 10, units = "cm")
rm(data_plot, plot)

Stark.

Head, K., Mayer, T. & Ries, J. (2010), “The erosion of colonial trade linkages after independence”. Journal of International Economics, 81(1):1-14 ↩︎
Es gibt mehrere Möglichkeiten diese section marks zu setzen, mehr dazu hier: https://support.rstudio.com/hc/en-us/articles/200484568-Code-Folding-and-Sections. ↩︎
Der Unterschied zu den naiven Regressionen der letzten Wochen ist, dass wir jetzt ja auch eine Zeitvariation haben. D.h. ↩︎
Andere Verteilungen sind auch möglich, so kann man Gravitationsmodelle auch bspw. mit einem Gamma Pseudo Maximum Likelihood Schätzer schätzen. Dafür dann family = Gamma(link = "log") setzen. ↩︎
Dieses Ergebnis tritt, interessanterweise regelmäßig auf, und ist ein “Puzzle” in der Forschung. ↩︎
Man muss nicht immer die Standardfarben oder Themes von ggplot2 nehmen.. ↩︎

Zuletzt aktualisiert am 19.05.2020