nunn_2008.R

library(haven)
library(ape)

# Load data
slave_trade_QJE <- read_dta("R/datasets/slave_trade_QJE.dta")

# Compute regression
mylm <- lm(ln_maddison_pcgdp2000~ln_export_area 
          # Colony fixed effects
          + colony0 + colony1 + colony2 + colony3 + colony4 + colony5 + colony6 + colony7 
          # Geographical controls
          + abs_latitude + longitude + rain_min + humid_max + low_temp + ln_coastline_area 
          + island_dum
          + legor_fr + region_n
          # + islam + legor_fr + legor_uk + region_n + region_s + region_e + region_w + region_c
          # Pre colonial prosperity controls
          + ln_avg_gold_pop + ln_avg_oil_pop + ln_avg_all_diamonds_pop
          + ln_pop_dens_1400
          ,data=slave_trade_QJE)

# Get some regression statistics
n <- nobs(mylm)
R2 <- summary(mylm)$adj.r.squared

# Compute standardized beta
r <- coef(summary(mylm))["ln_export_area", "Estimate"]
s <- coef(summary(mylm))["ln_export_area", "Std. Error"]
p <- coef(summary(mylm))["ln_export_area", 'Pr(>|t|)']

sigma_x <- sd(slave_trade_QJE$ln_export_area)
sigma_y <- sd(slave_trade_QJE$ln_maddison_pcgdp2000)

r_adjusted <- r * sigma_x / sigma_y
s_adjusted <- s * sigma_x / sigma_y

# Compute Moran's statistic
## Compute inverse distance matrix
slave_trade_QJE.dists <- as.matrix(dist(cbind(
  slave_trade_QJE$longitude, 
  slave_trade_QJE$abs_latitude)))

slave_trade_QJE.dists.inv <- 1/slave_trade_QJE.dists
diag(slave_trade_QJE.dists.inv) <- 0

## Compute Moran.I of residuals of GDP
moran_out <- Moran.I(resid(mylm), slave_trade_QJE.dists.inv)
moran_z <- (moran_out$observed - moran_out$expected) / moran_out$sd
moran_p <- moran_out$p.value

# Print results
print(sprintf('n = %d', n))
print(sprintf('r = %.2f (%.2f)', r, s))
print(sprintf('beta = %.2f (%.2f)', r_adjusted, s_adjusted))
print(sprintf('p = %.2f', p))
print(sprintf('Adjusted R2 = %.2f', R2))
print(sprintf("Moran's Z = %.2f", moran_z))
print(sprintf("Moran's p = %.2f", moran_p))