# 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

# SIR model of Covid-19 

# Parameters
# The user can change these parameters
# R0 <- 3.5      # Basic reproduction number (the average number of people
                 # infected by one infectious individual
# pgun <- 200    # When the policy is applied (days)
# petki <- 0     # Policy effectiveness
                 # 0 means full effect, 100 means no effect,
                 # After policy, R0 is reduced by petki percent, 
                 # i.e., R0 becomes R0*petki/100

# Constants
h_sure <- 15   # Duration of sickness, days
d_sure <- 15   # Duration of hospitalization, days
v_oran <- 10   # Death ratio, percent
N <- 1000      # Initial population
t <- 120       # İterasyon gün sayısı

# Cimulation function
covidSimul <- function(R0, pgun, petki) {
  # Set variables
  S <- I <- R <- V <- K <- rep(0, t)
  N <- rep(N, t)

  # Simulation starts...
  S[1] <- N[1] - 1
  I[1] <- 1
  
  for (i in c(2:t)) {
    pe <- ifelse(i >= pgun, petki / 100, 1) 
    ii <- ((R0 * pe)/h_sure)*S[i-1]*I[i-1]/(N[i-1] - V[i-1])
    rr <- I[i-1]/d_sure
    vv <- rr*(v_oran/100)
    kk <- rr - vv
    S[i] <- S[i-1] - ii
    I[i] <- I[i-1] + ii - rr
    R[i] <- R[i-1] + rr
    V[i] <- V[i-1] + vv
    K[i] <- K[i-1] + kk
  }
  dtab <- data.table(gun = c(1:t), S=S, I=I, R=R, V=V, K=K)
  return(dtab)
}

# Graphics function
simSekil <- function(veri, title, subtitle) {
  g1 <- ggplot(veri, aes(x = gun)) +
    geom_line(aes(y=S, color = "Healthy"), size = 1.5) +
    geom_line(aes(y=I, color = "Sick"), size = 1.5) +
    geom_line(aes(y=V, color = "Dead"), size = 1.5) +
    geom_line(aes(y=K, color = "Recovered"), size = 1.5) +
    labs(title = title, subtitle = subtitle,
         y = "Number of people", x = "Day") +
    theme_bw() +
    scale_color_manual(name = NULL,
                       values= c("Healthy" = "blue",
                                 "Sick" = "red",
                                 "Dead" = "black",
                                 "Recovered" = "orange"))
  return(g1)
}

# 
sim_results <- covidSimul(3.5, 200, 0)
simSekil(sim_results, "Covid-19 simulation results", "No policy")

sim_results <- covidSimul(3.5, 40, 30)
simSekil(sim_results, "Covid-19 simulation results", "Policy (30% effective) introduced on Day 40")