trade.R
In alpaca: Fit GLM's with High-Dimensional k-Way Fixed Effects

## ----eval = FALSE-------------------------------------------------------------
# # Required packages
# library(alpaca)
# library(data.table)
# library(haven)
# 
# # Import the data set
# cudata <- read_dta("dataaxj1.dta")
# setDT(cudata)
# 
# # Subset relevant variables
# var.nms <- c(
#   "exp1to2", "custrict11", "ldist", "comlang", "border", "regional",
#   "comcol", "curcol", "colony", "comctry", "cuwoemu", "emu", "cuc",
#   "cty1", "cty2", "year", "pairid"
#   )
# cudata <- cudata[, var.nms, with = FALSE]
# 
# # Generate identifiers required for structural gravity
# cudata[, pairid := factor(pairid)]
# cudata[, exp.time := interaction(cty1, year)]
# cudata[, imp.time := interaction(cty2, year)]
# 
# # Generate dummies for disaggregated currency unions
# cudata[, cuau := as.numeric(cuc == "au")]
# cudata[, cube := as.numeric(cuc == "be")]
# cudata[, cuca := as.numeric(cuc == "ca")]
# cudata[, cucf := as.numeric(cuc == "cf")]
# cudata[, cucp := as.numeric(cuc == "cp")]
# cudata[, cudk := as.numeric(cuc == "dk")]
# cudata[, cuea := as.numeric(cuc == "ea")]
# cudata[, cuec := as.numeric(cuc == "ec")]
# cudata[, cuem := as.numeric(cuc == "em")]
# cudata[, cufr := as.numeric(cuc == "fr")]
# cudata[, cugb := as.numeric(cuc == "gb")]
# cudata[, cuin := as.numeric(cuc == "in")]
# cudata[, cuma := as.numeric(cuc == "ma")]
# cudata[, cuml := as.numeric(cuc == "ml")]
# cudata[, cunc := as.numeric(cuc == "nc")]
# cudata[, cunz := as.numeric(cuc == "nz")]
# cudata[, cupk := as.numeric(cuc == "pk")]
# cudata[, cupt := as.numeric(cuc == "pt")]
# cudata[, cusa := as.numeric(cuc == "sa")]
# cudata[, cusp := as.numeric(cuc == "sp")]
# cudata[, cuua := as.numeric(cuc == "ua")]
# cudata[, cuus := as.numeric(cuc == "us")]
# cudata[, cuwa := as.numeric(cuc == "wa")]
# cudata[, cuwoo := custrict11]
# cudata[cuc %in% c("em", "au", "cf", "ec", "fr", "gb", "in", "us"), cuwoo := 0]
# 
# # Set missing trade flows to zero
# cudata[is.na(exp1to2), exp1to2 := 0]
# 
# # Re-scale trade flows
# cudata[, exp1to2 := exp1to2 / 1000]
# 
# # Construct binary and lagged dependent variable for the extensive margin
# cudata[, trade := as.numeric(exp1to2 > 0)]
# cudata[, ltrade := shift(trade), by = pairid]

## ----eval = FALSE-------------------------------------------------------------
# mod <- feglm(
#   exp1to2 ~ emu + cuwoo + cuau + cucf + cuec + cufr + cugb + cuin + cuus +
#     regional + curcol | exp.time + imp.time + pairid | cty1 + cty2 + year,
#   data   = cudata,
#   family = poisson()
#   )
# summary(mod, "sandwich")

## ----eval = FALSE-------------------------------------------------------------
# ## poisson - log link
# ##
# ## exp1to2 ~ emu + cuwoo + cuau + cucf + cuec + cufr + cugb + cuin +
# ##     cuus + regional + curcol | exp.time + imp.time + pairid |
# ##     cty1 + cty2 + year
# ##
# ## Estimates:
# ##           Estimate Std. error z value Pr(> |z|)
# ## emu       0.048895   0.010277   4.758  1.96e-06 ***
# ## cuwoo     0.765988   0.053272  14.379   < 2e-16 ***
# ## cuau      0.384469   0.118832   3.235   0.00121 **
# ## cucf     -0.125608   0.099674  -1.260   0.20760
# ## cuec     -0.877318   0.083451 -10.513   < 2e-16 ***
# ## cufr      2.095726   0.062952  33.291   < 2e-16 ***
# ## cugb      1.059957   0.034680  30.564   < 2e-16 ***
# ## cuin      0.169745   0.147029   1.154   0.24830
# ## cuus      0.018323   0.021530   0.851   0.39473
# ## regional  0.159181   0.008714  18.267   < 2e-16 ***
# ## curcol    0.386882   0.046827   8.262   < 2e-16 ***
# ## ---
# ## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# ##
# ## residual deviance= 35830.78,
# ## null deviance= 2245707.30,
# ## n= 1610165, l= [11227, 11277, 34104]
# ##
# ## ( 1363003 observation(s) deleted due to perfect classification )
# ##
# ## Number of Fisher Scoring Iterations: 20

## ----eval = FALSE-------------------------------------------------------------
# summary(mod, type = "clustered", cluster = ~ cty1 + cty2 + year)

## ----eval = FALSE-------------------------------------------------------------
# ## poisson - log link
# ##
# ## exp1to2 ~ emu + cuwoo + cuau + cucf + cuec + cufr + cugb + cuin +
# ##     cuus + regional + curcol | exp.time + imp.time + pairid |
# ##     cty1 + cty2 + year
# ##
# ## Estimates:
# ##          Estimate Std. error z value Pr(> |z|)
# ## emu       0.04890    0.09455   0.517   0.60507
# ## cuwoo     0.76599    0.24933   3.072   0.00213 **
# ## cuau      0.38447    0.22355   1.720   0.08546 .
# ## cucf     -0.12561    0.35221  -0.357   0.72137
# ## cuec     -0.87732    0.29493  -2.975   0.00293 **
# ## cufr      2.09573    0.30625   6.843  7.75e-12 ***
# ## cugb      1.05996    0.23766   4.460  8.19e-06 ***
# ## cuin      0.16974    0.30090   0.564   0.57267
# ## cuus      0.01832    0.05092   0.360   0.71898
# ## regional  0.15918    0.07588   2.098   0.03593 *
# ## curcol    0.38688    0.15509   2.495   0.01261 *
# ## ---
# ## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# ##
# ## residual deviance= 35830.78,
# ## null deviance= 2245707.30,
# ## n= 1610165, l= [11227, 11277, 34104]
# ##
# ## ( 1363003 observation(s) deleted due to perfect classification )
# ##
# ## Number of Fisher Scoring Iterations: 20

## ----eval = FALSE-------------------------------------------------------------
# library(car)
# h0_cus <- c("cuwoo", "cuau", "cucf", "cuec", "cufr", "cugb", "cuin", "cuus")
# linearHypothesis(
#   mod, h0_cus,
#   vcov. = vcov(mod, "clustered", cluster = ~ cty1 + cty2 + year)
#   )

## ----eval = FALSE-------------------------------------------------------------
# ## Linear hypothesis test
# ##
# ## Hypothesis:
# ## cuwoo = 0
# ## cuau = 0
# ## cucf = 0
# ## cuec = 0
# ## cufr = 0
# ## cugb = 0
# ## cuin = 0
# ## cuus = 0
# ##
# ## Model 1: restricted model
# ## Model 2: exp1to2 ~ emu + cuwoo + cuau + cucf + cuec + cufr + cugb + cuin +
# ##     cuus + regional + curcol | exp.time + imp.time + pairid |
# ##     cty1 + cty2 + year
# ##
# ## Note: Coefficient covariance matrix supplied.
# ##
# ##   Df  Chisq Pr(>Chisq)
# ## 1
# ## 2  8 96.772  < 2.2e-16 ***
# ## ---
# ## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

## ----eval = FALSE-------------------------------------------------------------
# mods <- feglm(
#   trade ~ cuwoemu + emu + regional | exp.time + imp.time + pairid,
#   data   = cudata,
#   family = binomial("logit")
#   )
# summary(mods)

## ----eval = FALSE-------------------------------------------------------------
# ## binomial - logit link
# ##
# ## trade ~ cuwoemu + emu + regional | exp.time + imp.time + pairid
# ##
# ## Estimates:
# ##          Estimate Std. error z value Pr(> |z|)
# ## cuwoemu   0.41037    0.05103   8.041  8.89e-16 ***
# ## emu       0.55108    0.23080   2.388   0.01696 *
# ## regional  0.10193    0.03271   3.116   0.00183 **
# ## ---
# ## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# ##
# ## residual deviance= 674579.11,
# ## null deviance= 1917254.65,
# ## n= 1384892, l= [11150, 11218, 29391]
# ##
# ## ( 1588276 observation(s) deleted due to perfect classification )
# ##
# ## Number of Fisher Scoring Iterations: 11

## ----eval = FALSE-------------------------------------------------------------
# modsbc <- biasCorr(mods, panel.structure = "network")
# summary(modsbc)

## ----eval = FALSE-------------------------------------------------------------
# ## binomial - logit link
# ##
# ## trade ~ cuwoemu + emu + regional | exp.time + imp.time + pairid
# ##
# ## Estimates:
# ##          Estimate Std. error z value Pr(> |z|)
# ## cuwoemu   0.37630    0.05107   7.369  1.72e-13 ***
# ## emu       0.47213    0.23056   2.048   0.04058 *
# ## regional  0.09942    0.03270   3.040   0.00236 **
# ## ---
# ## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# ##
# ## residual deviance= 674579.11,
# ## null deviance= 1917254.65,
# ## n= 1384892, l= [11150, 11218, 29391]
# ##
# ## ( 1588276 observation(s) deleted due to perfect classification )
# ##
# ## Number of Fisher Scoring Iterations: 11

## ----eval = FALSE-------------------------------------------------------------
# apesbc <- getAPEs(modsbc)
# summary(apesbc)

## ----eval = FALSE-------------------------------------------------------------
# ## Estimates:
# ##          Estimate Std. error z value Pr(> |z|)
# ## cuwoemu  0.015554   0.001957   7.946  1.92e-15 ***
# ## emu      0.019567   0.008212   2.383  0.017185 *
# ## regional 0.004079   0.001189   3.429  0.000605 ***
# ## ---
# ## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

## ----eval = FALSE-------------------------------------------------------------
# modd <- feglm(
#   trade ~ ltrade + cuwoemu + emu + regional | exp.time + imp.time + pairid,
#   data   = cudata,
#   family = binomial("logit")
#   )
# moddbc <- biasCorr(modd, L = 4, panel.structure = "network")
# summary(moddbc)

## ----eval = FALSE-------------------------------------------------------------
# ## binomial - logit link
# ##
# ## trade ~ ltrade + cuwoemu + emu + regional | exp.time + imp.time +
# ##     pairid
# ##
# ## Estimates:
# ##          Estimate Std. error z value Pr(> |z|)
# ## ltrade   2.166103   0.008003 270.663   < 2e-16 ***
# ## cuwoemu  0.268477   0.054482   4.928  8.31e-07 ***
# ## emu      0.312551   0.254091   1.230    0.2187
# ## regional 0.063077   0.035328   1.785    0.0742 .
# ## ---
# ## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# ##
# ## residual deviance= 602773.13,
# ## null deviance= 1899840.83,
# ## n= 1372471, l= [11054, 11116, 29299]
# ##
# ## ( 45048 observation(s) deleted due to missingness )
# ## ( 1555649 observation(s) deleted due to perfect classification )
# ##
# ## Number of Fisher Scoring Iterations: 11

## ----eval = FALSE-------------------------------------------------------------
# apedbc <- getAPEs(moddbc)
# summary(apedbc)

## ----eval = FALSE-------------------------------------------------------------
# ## Estimates:
# ##           Estimate Std. error z value Pr(> |z|)
# ## ltrade   0.1213377  0.0004681 259.231   < 2e-16 ***
# ## cuwoemu  0.0100223  0.0016486   6.079  1.21e-09 ***
# ## emu      0.0116944  0.0078456   1.491    0.1361
# ## regional 0.0023338  0.0010544   2.213    0.0269 *
# ## ---
# ## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1