pacman::p_load(readxl, tidyverse)
setwd("/Users/shanchaowang/Desktop/Submission Data/")

# read data
# GDP deflator from  1929 to 2020, with 2012 = 100
gdp = read_xls("BEA Implicit Price Index for GDP (2014-04-03).xls",range = "B6:CP8") %>%
  .[2,]
gdp = gather(gdp, key = "Year", value = "Deflator", `1929`:`2020`) %>% .[,-1]
# RD data in 2019 US$ rebase gdp to 2019 = 1
gdp$Deflator = gdp$Deflator/gdp$Deflator[gdp$Year == 2019]
gdp$Year = as.numeric(gdp$Year)

# value of production
VOP = read_xlsx("VA_State_US (2021-04-03).xlsx", range = "United States!B37:DJ37", col_names = F)
VOP = VOP[,-4] 
colnames(VOP) = 1910:2021
VOP = gather(VOP, key = "Year", value = "VOP") %>% as.data.frame()
VOP$Year = as.numeric(VOP$Year)

# merge VOP and gdp deflator
VOP = merge(VOP, gdp, by = "Year")
rm(gdp)

# deflate VOP in $1000, /1000 to convert to million $ in 2019 
VOP$VOP = VOP$VOP/VOP$Deflator/1000
VOP = VOP %>% select(Year, VOP) %>% rename(YEAR = Year)

# load constructed knowledge stock
load("Distribution.Rdata")
# load truncated distribution weights 
load("weights.RData")
weights_data = weights_data %>% rename(BLOOMEXPO = BloomExp,
                                       BLOOMREAL = BloomReal,
                                       BLOOMLINEAR = BloomLinear,
                                       BLOOMFLAT   = BloomFlat)
# read elasticity
e_OLS = read.csv("3_outputs/res_robust_ols_WI_std.csv") %>% select(model, elasticity, lower, upper) %>% filter(!is.na(lower)) 
e_CO  = read.csv("3_outputs/res_co_reg_WI_std.csv") %>% select(model, elasticity, lower, upper) %>% filter(!is.na(lower)) 
e_PW  = read.csv("3_outputs/res_pw_reg_WI_std.csv") %>% select(model, elasticity, lower, upper) %>% filter(!is.na(lower)) 
e_DOLS = read.csv("3_outputs/res_dols_sw_reg_WI_std.csv") %>% select(model, elasticity, lower, upper) %>% filter(!is.na(lower)) 
model.names = e_OLS$model

# BCR model
# Enter model name, discount rate, return estimated BCR and 95% confidence interval
mydata$BLOOMLINEAR_K = mydata$BLOOMLINEAR

BCR = function(model, r){
  BCR_data = merge(mydata %>% select('YEAR', model), VOP, by = "YEAR") %>% filter(YEAR >= 1956)
  Rt = BCR_data[which(BCR_data$YEAR == 1957), model] - BCR_data[which(BCR_data$YEAR == 1956), model]
  BCR_data = BCR_data %>% filter(YEAR >= 1957)
  if(model == "BLOOMLINEAR_K"){
    temp = BCR_data$VOP*sapply(BCR_data$YEAR-1957, function(x) (1+r)^(-x))/median(mydata$BLOOMLINEAR_K)
  }else{
    BCR_data$WEIGHT = weights_data[, model]
    temp = BCR_data$VOP/BCR_data[, model] * BCR_data$WEIGHT * sapply(BCR_data$YEAR-1957, function(x) (1+r)^(-x)) 
  }
  temp = sum(temp)
  out = cbind(model,
              e_OLS[which(e_OLS$model == model), 2:4] * temp,
              e_CO[which(e_CO$model == model), 2:4] * temp,
              e_PW[which(e_PW$model == model), 2:4] * temp,
              e_DOLS[which(e_DOLS$model == model), 2:4] * temp)
  colnames(out)[c(2, 5, 8, 11)] = c("OLS", "CO", "PW", "DOLS")
  return(out)
}

write.csv(lapply(model.names, BCR, r = 0.001) %>% bind_rows(), file = "4_outputs/BCR_001.csv", row.names = F)
write.csv(lapply(model.names, BCR, r = 0.03) %>% bind_rows(), file = "4_outputs/BCR_03.csv", row.names = F)
write.csv(lapply(model.names, BCR, r = 0.05) %>% bind_rows(), file = "4_outputs/BCR_05.csv", row.names = F)
write.csv(lapply(model.names, BCR, r = 0.1) %>% bind_rows(), file = "4_outputs/BCR_1.csv", row.names = F)