---
title: "esubt-analysis-comments"
output: html_document
date: "`r Sys.Date()`"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

if(!('devtools' %in% .packages(all.available = TRUE)))
  install.packages('devtools')
if(!('HARr' %in% .packages(all.available = TRUE)))
  devtools::install_git('https://github.com/USDA-ERS/MTED-HARr.git')
```


```{r include=FALSE }
# Read the past five versions of the GTAP database
  GTAPData = Map(function(f)HARr::read_har(sprintf('C:/Supplementary+materials/V7_flexAgg11Y%s/gsddat.har',f), FALSE, FALSE), c('04','07','11','14','17'))

# common script to calculate the tax on inputs and input/output values
  process = quote({
    tax = (VIFA + VDFA) / (VIFM + VDFM)
    output = apply(VIFA + VDFA, c(2, 3), sum) + apply(EVFA, c(2, 3), sum)
    outputMarket  = VDPM + VDGM + apply(VDFM, c(1,3),sum)+ apply(VXMD,c(1,2),sum)
    outputMarket[dimnames(VST)[[1]],] = outputMarket[dimnames(VST)[[1]],] + VST[,]
    outputTax = output
    outputTax[] = 1
    outputTax[dimnames(outputMarket)[[1]],] = outputMarket/output[dimnames(outputMarket)[[1]],]
    
    inputno = VIFM + VDFM
    input = VIFA + VDFA
    taxVa = apply(EVFA, c(2, 3), sum) / apply(VFM, c(2, 3), sum)
    va = apply(EVFA, c(2, 3), sum)
    vano = apply(VFM, c(2, 3), sum)
    regions = dimnames(VIFA)[[3]]
    commodities = dimnames(VIFA)[[2]]
    inputCommodities = dimnames(VIFA)[[1]]
  })
  
  GTAPDataProcessed = Map(function(f) within(f, eval(process)), GTAPData)
  

  combine = quote(do.call(rbind, Map(function(rr) {
    do.call(rbind, Map(function(cc) {
      rbind(
        do.call(rbind, Map(function(ii) {
          data.frame(
            inputType = 'input',
            inputName = ii,
            commodity = cc,
            region = rr,
            tax = tax[ii, cc, rr],
            outputTax = outputTax[cc, rr],
            input = input[ii, cc, rr],
            inputno = inputno[ii,cc,rr],
            output = output[cc, rr]
          )
        }, inputCommodities)),
        
        data.frame(
          inputType = 'factor',
          inputName = 'va',
          commodity = cc,
          region = rr,
          tax = taxVa[cc, rr],
          outputTax = outputTax[cc,rr],
          input = va[cc, rr],
          inputno = vano[cc,rr],
          output = output[cc, rr]
        )
        
      )
    }, commodities))
  }, regions)))
  
  cl = parallel::makeCluster(6)
  
  parallel::clusterExport(cl,'combine')
  
  GTAPDataCombined = parallel::clusterApplyLB(cl,GTAPDataProcessed, function(f)with(f,eval(combine)))
  
  parallel::stopCluster(cl)

  names(GTAPDataCombined) = names(GTAPDataProcessed)
  
  GTAPDataCombinedVersioned = Map(function(f,g){
    f$version = g
    return(f)
    }, GTAPDataCombined, order(names(GTAPDataCombined)))

    

  esubtData = do.call(rbind, GTAPDataCombinedVersioned)
  

  # remove all observation without noticeable taxation
  esubtDataClean =  esubtData[esubtData$tax > 1.0001 |  esubtData$tax < 1 / 1.0001, ]
  esubtDataClean =  esubtData[esubtData$tax > 1.0001 |  esubtData$tax < 1 / 1.0001, ]
  
  esubtDataClean = esubtDataClean[complete.cases(esubtDataClean),]
  
  # List of unique commodities (there are some extra ones because the coverage went up from 57 to 68)
 uc = unique(esubtDataClean$commodity)[!is.na(unique(esubtDataClean$commodity))]
  
  # Run the estimation for each commodity
  # Estimation is only based on the log of power of tax, log of value of output and fixed effects for commodities and regions
  
```

## IBN (Eq15)

```{r include=FALSE}
 # Selection of main crop
  uc <- c("pdr","wht","gro","v_f","osd","c_b","pfb","ocr")
  results_oneyear = Map(function(cc) {
      message(cc)
      mm = lm(
        formula = I(log(input) - log(output)) ~   I(log(tax)+log(outputTax)) +  inputName + region + version,
        data = esubtDataClean[esubtDataClean$commodity == cc, ]
      )
    }, uc)
    
  names(results_oneyear) = uc
    
  resultsSummaries = Map(function(f)
    list(coefficients = summary(f)$coefficients, r2 = summary(f)$adj.r.squared,df=summary(f)$df[2], k=summary(f)$df[1],
    lf= (summary(f)$df[2]-1)*((summary(f)$df[1]+summary(f)$df[2])^(1/(summary(f)$df[1]+summary(f)$df[2]))-1)), results_oneyear)
    
  results = resultsSummaries
  table = do.call(cbind,Map(function(mm,com){
    c(Estimate = round(1-mm$coefficients['I(log(tax) + log(outputTax))',1],3),
    Standard = round(mm$coefficients['I(log(tax) + log(outputTax))',2],3),
    Significanced   = 2*round(pnorm(-abs(1 - mm$coefficients['I(log(tax) + log(outputTax))',1]) / mm$coefficients['I(log(tax) + log(outputTax))',2]),3))}, results[sort(names(results))], sort(names(results))))
  tableIBN <- matrix(c(table),ncol = 1)
  colnames(tableIBN) <- "IBN (Eq15)"
  table
```


## CRS with all inputs (Eq18)
```{r}
  #######################################################
  ## tackling CRS issue on all inputs: Second approach log(vifm+vdfm) = -sigma*log(tax)+dummies
# Selection of main crop
  uc <- c("pdr","wht","gro","v_f","osd","c_b","pfb","ocr")
  results_oneyear = Map(function(cc) {
      message(cc)
      mm = lm(
        formula = I(log(inputno))  ~   I(log(tax))  +  inputName + region + version,
        data = esubtDataClean[esubtDataClean$commodity == cc, ]
      )
    }, uc)
    
  names(results_oneyear) = uc
    
  resultsSummaries = Map(function(f)
    list(coefficients = summary(f)$coefficients, r2 = summary(f)$adj.r.squared,df=summary(f)$df[2], k=summary(f)$df[1],
    lf= (summary(f)$df[2]-1)*((summary(f)$df[1]+summary(f)$df[2])^(1/(summary(f)$df[1]+summary(f)$df[2]))-1)), results_oneyear)
    
  results = resultsSummaries
  table = do.call(cbind,Map(function(mm,com){
    c(Estimate = round(-mm$coefficients['I(log(tax))',1],3),
    Standard = round(mm$coefficients['I(log(tax))',2],3),
    Significance   = 2*round(pnorm(- abs( - mm$coefficients['I(log(tax))',1])/ mm$coefficients['I(log(tax))',2]),3))}, results[sort(names(results))], sort(names(results))))
  tableCRS <- matrix(c(table),ncol = 1)
  colnames(tableCRS) <- "CRS (Eq18)"
  table
```
## IBN with Key Inputs (Eq15)

```{r}
  #####################################
  # Selection of main crop
  uc <- c("pdr","wht","gro","v_f","osd","c_b","pfb","ocr")
  results_oneyear = Map(function(cc) {
      message(cc)
      ic <- c("chm","ely","p_c","va", cc) # selection of inputs
      mm = lm(
        formula = I(log(input) - log(output)) ~   I(log(tax)+log(outputTax)) +  inputName + region + version,
        data = esubtDataClean[esubtDataClean$commodity == cc & esubtDataClean$inputName %in% ic , ]
      )
    }, uc)
    
  names(results_oneyear) = uc
    
  resultsSummaries = Map(function(f)
    list(coefficients = summary(f)$coefficients, r2 = summary(f)$adj.r.squared,df=summary(f)$df[2], k=summary(f)$df[1],
    lf= summary(f)$df[2]*((summary(f)$df[1]+summary(f)$df[2])^(1/(summary(f)$df[1]+summary(f)$df[2]))-1)), results_oneyear)
  
  results = resultsSummaries
  table = do.call(cbind,Map(function(mm,com){
    c(Estimate = round(1-mm$coefficients['I(log(tax) + log(outputTax))',1],3),
    Standard = round(mm$coefficients['I(log(tax) + log(outputTax))',2],3),
    Significance   = 2*round(pnorm(- abs(1 - mm$coefficients['I(log(tax) + log(outputTax))',1]) / mm$coefficients['I(log(tax) + log(outputTax))',2]),3))}, results[sort(names(results))], sort(names(results))))
  tableIBN_INP <- matrix(c(table),ncol = 1)
  colnames(tableIBN_INP) <- "Key inputs (Eq15)"
  table
```


## CRS on Key Inputs (Eq18)
```{r}
  ############################################################
# Selection of main crop
  uc <- c("pdr","wht","gro","v_f","osd","c_b","pfb","ocr")
  results_oneyear = Map(function(cc) {
      message(cc)
    ic <- c("chm","ely","p_c","va", cc) # selection of inputs
      mm = lm(
        formula = I(log(inputno))  ~   I(log(tax))  +  inputName + region + version,
        data = esubtDataClean[esubtDataClean$commodity == cc & esubtDataClean$inputName %in% ic , ]
      )
    }, uc)
    
  names(results_oneyear) = uc
    
  resultsSummaries = Map(function(f)
    list(coefficients = summary(f)$coefficients, r2 = summary(f)$adj.r.squared,df=summary(f)$df[2], k=summary(f)$df[1],
    lf= (summary(f)$df[2]-1)*((summary(f)$df[1]+summary(f)$df[2])^(1/(summary(f)$df[1]+summary(f)$df[2]))-1)), results_oneyear)
    
  results = resultsSummaries
  table = do.call(cbind,Map(function(mm,com){
    c(Estimate = round(-mm$coefficients['I(log(tax))',1],3),
    Standard = round(mm$coefficients['I(log(tax))',2],3),
    Significance   = 2*round(pnorm( -abs(- mm$coefficients['I(log(tax))',1])/ mm$coefficients['I(log(tax))',2]),3))}, results[sort(names(results))], sort(names(results))))
  tableCRS_INP <- matrix(c(table),ncol = 1)
  colnames(tableCRS_INP) <- "Key inputs (Eq18)"
  table
```
  
## Table 1
```{r}
  ######################
  name1 <- matrix(c(rep("Sugar crops",3),rep("Coarse grains",3),rep("Other crops",3),
                    rep("Oilseeds",3),rep("Paddy rice",3),rep("Plant-based fibers",3),
                    rep("Vegetable and fruits",3),rep("Wheat",3)),ncol = 1)
  name2 <- matrix(rep(c("Estimate","Standard error", "Significance"),8),ncol = 1)
  
  table1 <- data.frame(name1, name2, tableIBN, tableCRS, tableIBN_INP, tableCRS_INP)
  colnames(table1) <- c("", "", "IBN (Eq15)", "CRS (Eq18)", "Key inputs (Eq15)", "Key inputs (Eq18)")
  table1
```
 
 
# Dynamic

```{r}
# Calculate the long-run elasticities
  esubtData$previousVersion = esubtData$version - 1

  esubtDataMerged = merge(esubtData,esubtData, by.x=c('inputType','inputName','commodity','region','version'), by.y=c('inputType','inputName','commodity','region','previousVersion'))

  esubtDataMergedClean = esubtDataMerged[esubtDataMerged$tax.y>1.0001 | esubtDataMerged$tax.y<1/1.0001,]
  esubtDataMergedClean=esubtDataMergedClean[complete.cases(esubtDataMergedClean),]

```

## IBN (Eq19)

```{r}
  uc <- c("pdr","wht","gro","v_f","osd","c_b","pfb","ocr")
  
  results_threeyears = Map(function(cc) {
    message(cc)
    mm = lm(I(log(input.y)-log(output.y))~ + log(tax.y/tax.x * outputTax.y/outputTax.x) +region + inputName + version, data = esubtDataMergedClean[esubtDataMergedClean$commodity==cc,])
  }, uc)
  
  names(results_threeyears) = uc
  
  
  resultsSummariesLR = Map(function(f)
      list(coefficients = summary(f)$coefficients, r2 = summary(f)$adj.r.squared,df=summary(f)$df[2], k=summary(f)$df[1],
      lf= (summary(f)$df[2]-1)*((summary(f)$df[1]+summary(f)$df[2])^(1/(summary(f)$df[1]+summary(f)$df[2]))-1)), results_threeyears)
    
  results = resultsSummariesLR
  table = do.call(cbind,Map(function(mm,com){
    c(Estimate = round(1-mm$coefficients['log(tax.y/tax.x * outputTax.y/outputTax.x)',1],3),
    Standard = round(mm$coefficients['log(tax.y/tax.x * outputTax.y/outputTax.x)',2],3),
    Significance   = 2* round(pnorm(- abs(1 - mm$coefficients['log(tax.y/tax.x * outputTax.y/outputTax.x)',1]) / mm$coefficients['log(tax.y/tax.x * outputTax.y/outputTax.x)',2]),3))}, results[sort(names(results))], sort(names(results))))
  table_D_IBN <- matrix(c(table),ncol = 1)
  colnames(table_D_IBN) <- "IBN (Eq19)"
  table
```

## Our specifcation with all inputs (Eq25)
```{r}
  uc <- c("pdr","wht","gro","v_f","osd","c_b","pfb","ocr")
  
  results_threeyears = Map(function(cc) {
    message(cc)
    mm = lm(I(log(input.y)-log(output.y))~ + log(tax.y/tax.x * outputTax.y/outputTax.x) + log(tax.x * outputTax.x) + 
              region + inputName + version, data = esubtDataMergedClean[esubtDataMergedClean$commodity==cc,])
  }, uc)

  names(results_threeyears) = uc
  
  resultsSummariesLR = Map(function(f)
      list(coefficients = summary(f)$coefficients, r2 = summary(f)$adj.r.squared,df=summary(f)$df[2], k=summary(f)$df[1],
      lf= (summary(f)$df[2]-1)*((summary(f)$df[1]+summary(f)$df[2])^(1/(summary(f)$df[1]+summary(f)$df[2]))-1)),
       results_threeyears)
    
  results = resultsSummariesLR
  table = do.call(cbind,Map(function(mm,com){
    c(Estimate = round(1-mm$coefficients['log(tax.x * outputTax.x)',1],3),
    Standard = round(mm$coefficients['log(tax.x * outputTax.x)',2],3),
    Significance   = 2 * round(pnorm(-abs(1 - mm$coefficients['log(tax.x * outputTax.x)',1]) / mm$coefficients['log(tax.x * outputTax.x)',2]),3))}, results[sort(names(results))], sort(names(results))))
  table_D_IBN_SP <- matrix(c(table),ncol = 1)
  colnames(table_D_IBN_SP) <- "Our specification (Eq25)"
  table

```


## Our specification on relevent inputs (Eq25)
```{r}

  uc <- c("pdr","wht","gro","v_f","osd","c_b","pfb","ocr")
   results_threeyears = Map(function(cc) {
    message(cc)
     ic <- c("chm","ely","p_c","va", cc) # selection of inputs
    mm = lm(I(log(input.y)-log(output.y))~ + log(tax.y/tax.x * outputTax.y/outputTax.x) + log(tax.x * outputTax.x) + 
              region + inputName + version, data = esubtDataMergedClean[esubtDataMergedClean$commodity==cc & esubtDataMergedClean$inputName %in% ic,])
  
  }, uc)

  names(results_threeyears) = uc
  
  resultsSummariesLR = Map(function(f)
      list(coefficients = summary(f)$coefficients, r2 = summary(f)$adj.r.squared,df=summary(f)$df[2], k=summary(f)$df[1],
      lf= (summary(f)$df[2]-1)*((summary(f)$df[1]+summary(f)$df[2])^(1/(summary(f)$df[1]+summary(f)$df[2]))-1)), results_threeyears)
    
  results = resultsSummariesLR
  table = do.call(cbind,Map(function(mm,com){
    c(Estimate = round(1-mm$coefficients['log(tax.x * outputTax.x)',1],3),
    Standard = round(mm$coefficients['log(tax.x * outputTax.x)',2],3),
    Significance   = 2 * round(pnorm( - abs(1 - mm$coefficients['log(tax.x * outputTax.x)',1]) / mm$coefficients['log(tax.x * outputTax.x)',2]),3))}, results[sort(names(results))], sort(names(results))))
  table_D_IBN_SP_INP <- matrix(c(table),ncol = 1)
  colnames(table_D_IBN_SP_INP) <- "Key inputs (Eq25)"
  table
```

## Table2
```{r}
  name1 <- matrix(c(rep("Sugar crops",3),rep("Coarse grains",3),rep("Other crops",3),
                    rep("Oilseeds",3),rep("Paddy rice",3),rep("Plant-based fibers",3),
                    rep("Vegetable and fruits",3),rep("Wheat",3)),ncol = 1)
  name2 <- matrix(rep(c("Estimate","Standard error", "Significance"),8),ncol = 1)
  
  table2 <- data.frame(name1,name2,table_D_IBN,table_D_IBN_SP,table_D_IBN_SP_INP)
  colnames(table2) <- c("", "", "IBN (Eq19)", "Our specification (Eq25)", "Key inputs (Eq25)")
  table2
```