## packages
library("lmtest")
library("sandwich")

## read data
DJ <- read.table("data.dj", header = TRUE, na.strings = c("-999.0", "-999.00"))


############################
## Durlauf & Johnson (DJ) ##
############################

## use only non-oil countries
dj <- subset(DJ, NONOIL == 1)
## rescale only population growth
dj$POPGRO <- dj$POPGRO/100

## model equation (taken from MRW, Table V)
dj_model <-  I(log(GDP85) - log(GDP60)) ~ log(GDP60) + log(IONY) + log(POPGRO + 0.05) + log(SCHOOL)
## but MRW used different scaling
mrw_model <- I(log(GDP85) - log(GDP60)) ~ log(GDP60) + log(IONY/100) + log(POPGRO + 0.05) + log(SCHOOL/100)

## Table II
dj_mrw <- lm(mrw_model, data = dj)
dj_sub1 <- lm(dj_model, data = dj, subset = GDP60 < 1800 & LIT60 < 50)
dj_sub2 <- lm(dj_model, data = dj, subset = GDP60 >= 1800 & LIT60 >= 50)
coeftest(dj_mrw)
coeftest(dj_sub1, vcov = sandwich)
coeftest(dj_sub2, vcov = sandwich)

## Table V
dj_node1 <- lm(dj_model, data = dj, subset = GDP60 < 800)
dj_node2 <- lm(dj_model, data = dj, subset = GDP60 >= 800 & LIT60 < 46)
dj_node3 <- lm(dj_model, data = dj, subset = GDP60 >= 800 & LIT60 >= 46 & GDP60 < 4850)
dj_node4 <- lm(dj_model, data = dj, subset = GDP60 >= 800 & LIT60 >= 46 & GDP60 >= 4850)
coeftest(dj_node1)
coeftest(dj_node2, vcov = sandwich)
coeftest(dj_node3, vcov = sandwich)
coeftest(dj_node4)

## summary: reproducible (with some effort)
##   - scaling of regressors is not obvious
##     (and different for full-sample MRW model)
##   - type of standard errors inconsistent
##   - wrong thresholds in Table II header


###############################
## Mankiw, Romer, Weil (MRW) ##
###############################

## use full data set
## rescale all variables given in percent
mrw <- DJ
mrw$IONY <- mrw$IONY/100
mrw$POPGRO <- mrw$POPGRO/100
mrw$SCHOOL <- mrw$SCHOOL/100

## Table I -> see MP (mp_cd1)
mrw_1n <- lm(log(GDP85) ~ log(IONY) + log(POPGRO + 0.05), data = mrw, subset = NONOIL == 1)
mrw_1i <- lm(log(GDP85) ~ log(IONY) + log(POPGRO + 0.05), data = mrw, subset = INTER == 1)
mrw_1o <- lm(log(GDP85) ~ log(IONY) + log(POPGRO + 0.05), data = mrw, subset = OECD == 1)
summary(mrw_1n)
summary(mrw_1i)
summary(mrw_1o)

## Table II -> see MP (mp_cd2)
mrw_2n <- lm(log(GDP85) ~ log(IONY) + log(POPGRO + 0.05) + log(SCHOOL), data = mrw, subset = NONOIL == 1)
mrw_2i <- lm(log(GDP85) ~ log(IONY) + log(POPGRO + 0.05) + log(SCHOOL), data = mrw, subset = INTER == 1)
mrw_2o <- lm(log(GDP85) ~ log(IONY) + log(POPGRO + 0.05) + log(SCHOOL), data = mrw, subset = OECD == 1)
summary(mrw_2n)
summary(mrw_2i)
summary(mrw_2o)

## Table III
mrw_3n <- lm(I(log(GDP85) - log(GDP60)) ~ log(GDP60), data = mrw, subset = NONOIL == 1)
mrw_3i <- lm(I(log(GDP85) - log(GDP60)) ~ log(GDP60), data = mrw, subset = INTER == 1)
mrw_3o <- lm(I(log(GDP85) - log(GDP60)) ~ log(GDP60), data = mrw, subset = OECD == 1)
summary(mrw_3n)
summary(mrw_3i)
summary(mrw_3o)

## Table IV
mrw_4n <- lm(I(log(GDP85) - log(GDP60)) ~ log(GDP60) + log(IONY) + log(POPGRO + 0.05),
  data = mrw, subset = NONOIL == 1)
mrw_4i <- lm(I(log(GDP85) - log(GDP60)) ~ log(GDP60) + log(IONY) + log(POPGRO + 0.05),
  data = mrw, subset = INTER == 1)
mrw_4o <- lm(I(log(GDP85) - log(GDP60)) ~ log(GDP60) + log(IONY) + log(POPGRO + 0.05),
  data = mrw, subset = OECD == 1)
summary(mrw_4n)
summary(mrw_4i)
summary(mrw_4o)

## Table V -> see DJ (dj_mrw)
mrw_5n <- lm(I(log(GDP85) - log(GDP60)) ~ log(GDP60) + log(IONY) + log(POPGRO + 0.05) + log(SCHOOL),
  data = mrw, subset = NONOIL == 1)
mrw_5i <- lm(I(log(GDP85) - log(GDP60)) ~ log(GDP60) + log(IONY) + log(POPGRO + 0.05) + log(SCHOOL),
  data = mrw, subset = INTER == 1)
mrw_5o <- lm(I(log(GDP85) - log(GDP60)) ~ log(GDP60) + log(IONY) + log(POPGRO + 0.05) + log(SCHOOL),
  data = mrw, subset = OECD == 1)
summary(mrw_5n)
summary(mrw_5i)
summary(mrw_5o)

## Table VI
mrw_6n <- lm(I(log(GDP85) - log(GDP60)) ~ log(GDP60) + I(log(IONY) - log(POPGRO + 0.05)) +
  I(log(SCHOOL) - log(POPGRO + 0.05)), data = mrw, subset = NONOIL == 1)
mrw_6i <- lm(I(log(GDP85) - log(GDP60)) ~ log(GDP60) + I(log(IONY) - log(POPGRO + 0.05)) +
  I(log(SCHOOL) - log(POPGRO + 0.05)), data = mrw, subset = INTER == 1)
mrw_6o <- lm(I(log(GDP85) - log(GDP60)) ~ log(GDP60) + I(log(IONY) - log(POPGRO + 0.05)) +
  I(log(SCHOOL) - log(POPGRO + 0.05)), data = mrw, subset = OECD == 1)
summary(mrw_6n)
summary(mrw_6i)
summary(mrw_6o)

## summary: reproducible (provided regressors were scaled correctly)
##   - small deviations for many coefficients and standard errors


###################################
## Masanjala & Papageorgiou (MP) ##
###################################

## use only non-oil countries
mp <- subset(DJ, NONOIL == 1)
## use unscaled data directly

## Table I (based on PWT 4.0)
mp_cd1 <- lm(log(GDP85) ~ log(IONY) + log(POPGRO + 5), data = mp)
mp_ces1 <- lm(log(GDP85) ~ log(IONY) + log(POPGRO + 5) + I((log(IONY) - log(POPGRO + 5))^2), data = mp)
mp_cd2 <- lm(log(GDP85) ~ log(IONY) + log(POPGRO + 5) + log(SCHOOL), data = mp)
mp_ces2 <- lm(log(GDP85) ~ log(IONY) + log(POPGRO + 5) + log(SCHOOL) +
  I((log(IONY) - log(POPGRO + 5))^2) + I((log(SCHOOL) - log(POPGRO + 5))^2) +
  I((log(IONY) - log(SCHOOL))^2), data = mp)
coeftest(mp_cd1, vcov = sandwich)
coeftest(mp_ces1, vcov = sandwich)
coeftest(mp_cd2, vcov = sandwich)
coeftest(mp_ces2, vcov = sandwich)

## summary: reproducible (no scaling of regressors required)