# Load libraries
library(data.table)
library(haven)
library(ggplot2)
library(gganimate)

# Download the data
pwt <- read_stata("http://www.rug.nl/ggdc/docs/pwt90.dta")

# Make the data a data.table object
pwt <- data.table(pwt)

# Select a subset of variables
pwts <- pwt[, list(countrycode, year, rgdpo, pop, emp, avh, hc, rgdpna, rkna)]

# Structure
str(pwts)

# Generate GDP per capita variable - 000 USD
pwts[, gdpc := log(rgdpo / pop)]

# Drop the outlier
pwts <- pwts[!countrycode == "BMU"]

ggplot(pwts, aes(x=hc, y=gdpc, size=pop)) + geom_point()

pwts[, tr := countrycode == "TUR"]

ggplot(pwts, aes(x=hc, y=gdpc, size=pop, color = tr)) + geom_point()

ggplot(pwts[year %in% seq(from=1950, to=2010, by = 5)], 
       aes(x=hc, y=gdpc, size=pop, color = tr)) +
  geom_point() + facet_wrap(~year, ncol=3)

g1 <- ggplot(pwts, aes(x = hc, y=gdpc, size=pop, color = tr)) + 
  geom_point() + theme(legend.position="none") +
  labs(title = "Year: {round(frame_time, 0)}", x = "Education level", 
       y = "GDP per capita (log)") +
  transition_time(year)

g1

T <- 100    # Simulation period
alfa <- 0.3 # Output elasticity of capital
s <- 0.20   # Saving rate
n <- 0.02   # Population growth rate
d <- 0.05   # Depreciation rate
t <- 0.02   # Rate of technological change
ts <- 0.02   # Standard deviation of technological change

set.seed(123)
dat <- data.table(t = c(1:T), K = 1, L = 1, A = 1)
dat[t == 1, Q := A*(K^alfa)*(L^(1-alfa))]

# Iterations
for (i in c(2:T)) {
  dat$L[i] <- dat$L[i-1]*(1+n)
  dat$K[i] <- (dat$K[i-1] * (1-d)) + s * dat$Q[i-1]
  dat$A[i] <- dat$A[i-1] * (1 + rnorm(1, t, ts))
  dat$Q[i] <- (dat$K[i]^alfa) * (dat$A[i] * dat$L[i])^(1-alfa)
}

getGR <- function(x) x/shift(x) - 1

ggplot(dat, aes(x = K/(L*A), y = Q/(L*A))) + geom_point() + 
  geom_hline(yintercept = 0) + 
  labs(title = "Labor productivity and capital intensity",
       x = "Capital intensity (in efficiency units)",
       y = "Labor productivity (in efficiency units)")

ggplot(dat, aes(x = t, y = getGR(A))) + 
  geom_hline(yintercept = 0) + geom_line() + 
  labs(title = "Productivity growth rate",
       x = "Time", y = "Productivity growth rate (%)")

ggplot(dat, aes(x = t, y = getGR(Q))) + geom_line() + 
  geom_hline(yintercept = 0) + 
  labs(title = "Output growth rate",
       x = "Time", y = "Output growth rate (%)")

ggplot(dat, aes(x = t, y = getGR(Q/L))) + geom_line() + 
  geom_hline(yintercept = 0) + 
  labs(title = "Labor productivity growth rate",
       x = "Time", y = "Labor productivity growth rate (%)")

ggplot(dat, aes(x = t, y = getGR(K/L))) + 
  geom_hline(yintercept = 0) + geom_line() + 
  labs(title = "Capital intensity growth rate",
       x = "Time", y = "Capital intensity growth rate (%)")

solow <- function(model = "m0", T = 100, alfa = 0.3, s = 0.2, n = 0.02, d = 0.05, t = 0.02, ts = 0.02) {

  # T <- 100    # Simulation period
  # alfa <- 0.3 # Output elasticity of capital
  # s <- 0.20   # Saving rate
  # n <- 0.02   # Population growth rate
  # d <- 0.05   # Depreciation rate
  # t <- 0.02   # Rate of technological change
  # ts <- 0.02   # Standard deviation of technological change

  set.seed(123)
  dat <- data.table(t = c(1:T), K = 1, L = 1, A = 1, model = model)
  dat[t == 1, Q := A*(K^alfa)*(L^(1-alfa))]

  # Iterations
  for (i in c(2:T)) {
    dat$L[i] <- dat$L[i-1]*(1+n)
    dat$K[i] <- (dat$K[i-1] * (1-d)) + s * dat$Q[i-1]
    dat$A[i] <- dat$A[i-1] * (1 + rnorm(1, t, ts))
    dat$Q[i] <- (dat$K[i]^alfa) * (dat$A[i] * dat$L[i])^(1-alfa)
  }
  return(dat)
}

m1 <- solow(t = 0,    ts = 0,   model = "m1")
m2 <- solow(t = 0.02, ts = 0,    model = "m2")
m3 <- solow(t = 0.02, ts = 0.02, model = "m3")

models <- rbindlist(list(m1, m2, m3))

getGR <- function(x) {
  g <- x/shift(x) - 1
  g[dat$t == 1] <- NA
  return(g)
}

ggplot(models, aes(x = t, y = getGR(Q), color = model)) + 
  geom_line(size = 1.1) + theme_bw() +
  geom_hline(yintercept = 0) + 
  labs(title = "Output growth rate",
       x = "Time", y = "Output growth rate (%)")

ggplot(models, aes(x = t, y = getGR(Q/L), color = model)) + 
  geom_line(size = 1.1) + theme_bw() +
  geom_hline(yintercept = 0) + 
  labs(title = "Labor productivity growth rate",
       x = "Time", y = "Labor productivity growth rate (%)")

ggplot(models, aes(x = t, y = getGR(K/L), color = model)) + 
  geom_line(size = 1.1) + theme_bw() +
  geom_hline(yintercept = 0) + geom_line() + 
  labs(title = "Capital intensity growth rate",
       x = "Time", y = "Capital intensity growth rate (%)")