library("foreign")
library("sandwich")
library("lmtest")
install.packages("KRLS")
library("KRLS")
growth <- read.dta("growthdata.dta")
growth <- growth[order(growth$yearsschool), ]
### Basic Fits
## OLS
l.out <- lm(growth ~ yearsschool, data = growth)
coeftest(l.out, vcovHC(l.out, type = "HC1"))
## KRLS
library("KRLS")
growth <- read.dta("growthdata.dta")
growth <- growth[order(growth$yearsschool), ]
k.out <- krls(y = growth$growth, X = growth$yearsschool)
summary(k.out)
## OLS 3rd order poly
l.out2 <- lm(growth ~ yearsschool + I(yearsschool^2) + I(yearsschool^3),
data = growth)
## Figure 1: Fitted values from KRLS and OLS models
plot(y = growth$growth, x = growth$yearsschool)
lines(y = l.out$fitted, x = growth$yearsschool, col = "green")
lines(y = k.out$fitted, x = growth$yearsschool, col = "red")
lines(y = l.out2$fitted, x = growth$yearsschool, col = "orange", lty = 2)
legend("bottomright",
legend = c("KRLS fit", "OLS fit", "OLS 3rd Poly Fit"),
col = c("red", "green", "orange"),
lty = c(1, 1, 2))
### Pointwise Partial Derivatives
## OLS
l.out <- lm(growth ~ yearsschool + assassinations,
data = growth)
coeftest(l.out, vcovHC(l.out, type = "HC1"))
## KRLS
k.out <- krls(y = growth$growth,
X = growth[c("yearsschool", "assassinations")])
summary(k.out)
## Figure 2: Conditional relationship between growth and schooling (controlling for assassinations).
CRu <- l.out$resid + growth$yearsschool * coef(l.out)[2]
plot(CRu, x = growth$yearsschool,
ylab = "Component + Residual:Growth|Assassinations",
xlab = "average years of Schooling")
lines(lowess(CRu ~ growth$yearsschool))
abline(lm(CRu ~ growth$yearsschool))
## Heterogeneity of the marginal effect in KRLS
summary(k.out$derivatives[, "yearsschool"])
## Figure 3: Distribution of pointwise marginal effect of schooling on growth
hist(k.out$derivatives[, "yearsschool"], xlab = "d_yearsschool", freq = FALSE, main = NA)
## Figure 4: Pointwise marginal effect of schooling and level of schooling
plot(y = k.out$derivatives[, "yearsschool"],
x = growth$yearsschool)
lines(lowess(k.out$derivatives[, "yearsschool"] ~ growth$yearsschool))
### The Full Model
l.out <- lm(growth ~ rgdp60 + tradeshare + yearsschool + assassinations,
data = growth)
coeftest(l.out, vcovHC(l.out, type = "HC1"))
covars <- c("rgdp60", "tradeshare",
"yearsschool", "assassinations")
k.out <- krls(y = growth$growth,
X = growth[, covars])
summary(k.out)
## fixing the schooling
l.out2 <- lm(growth ~ rgdp60 + tradeshare +
yearsschool + I(yearsschool^2) + I(yearsschool^3) +
assassinations,
data = growth)
coeftest(l.out2, vcovHC(l.out2, type = "HC1"))
## heterogeneity in trade effect
growth$d_tradeshare <- k.out$derivatives[, "tradeshare"]
summary(lm(d_tradeshare ~ rgdp60 + tradeshare + yearsschool + assassinations,
data = growth))
## Figure 5: Pointwise marginal effect of trade share and level of trade share.
plot(growth$d_tradeshare ~ growth$tradeshare, xlab = "trade share of GDP")
lines(lowess(growth$d_tradeshare ~ growth$tradeshare))
## Figure 6: Pointwise marginal effect of trade share and number of assassinations.
plot(growth$d_tradeshare ~ growth$assassinations, xlab = "number of assassinations")
lines(lowess(growth$d_tradeshare ~ growth$assassinations))
## fixing the other predictors
l.out3 <- lm(growth ~ rgdp60 + tradeshare + I(tradeshare^2) + I(tradeshare^3) +
assassinations + I(tradeshare*assassinations) +
I(tradeshare^2*assassinations) + I(tradeshare^3*assassinations) +
yearsschool + I(yearsschool^2) + I(yearsschool^3),
data = growth)
coeftest(l.out3, vcovHC(l.out3, type = "HC1"))
## augmented ols fits almost as well as krls.
cor(cbind(k.out$fitted, l.out3$fitted, l.out$fitted))
### LOO error of krls, lm, and the augmented lm
covars <- c("rgdp60", "tradeshare", "yearsschool", "assassinations")
LOOerrors <- matrix(NA, nrow = nrow(growth), ncol = 3)
colnames(LOOerrors) <- c("krls", "ols1", "ols2")
for (leftout in 1:nrow(growth)) {
krls.out <- krls(y = growth$growth[-leftout], X = growth[-leftout, covars], derivative = FALSE)
LOOerrors[leftout, 1] <- growth$growth[leftout] - predict(krls.out, newdata = growth[leftout, covars])$fit
lm.out1 <- lm(growth ~ rgdp60 + tradeshare + assassinations + yearsschool, data = growth[-leftout, ])
LOOerrors[leftout, 2] <- growth$growth[leftout]-predict(lm.out1, newdata = growth[leftout, ])
lm.out2 <- lm(growth ~ rgdp60 + tradeshare + I(tradeshare^2) + I(tradeshare^3) + assassinations +
I(tradeshare*assassinations) + I(tradeshare^2*assassinations) + I(tradeshare^3*assassinations) +
yearsschool + I(yearsschool^2) + I(yearsschool^3), data = growth[-leftout, ])
LOOerrors[leftout, 3] <- growth$growth[leftout] - predict(lm.out2, newdata = growth[leftout, ])
}
## Get MSEs
apply(LOOerrors, 2, function(x) mean(x^2))
### Binary Predictors
growth$yearsschool3 <- growth$yearsschool>3
k.out <- krls(y = growth$growth,
X = growth[, c("rgdp60", "tradeshare", "yearsschool3", "assassinations")])
summary(k.out)