R/StateiorFunctions.R

In USEPA/useeior: USEEIO R modeling software

# Functions for handling data from stateior, https://github.com/USEPA/stateior

#' Load two-region IO data of model iolevel and year from user's local directory
#' or the EPA Data Commons.
#' @description Load two-region IO data of model iolevel and year from user's
#' local directory or the EPA Data Commons.
#' @param model An EEIO form USEEIO model object with model specs and IO meta data loaded.
#' @param dataname Name of desired IO data, can be "Make", "Use", "DomesticUse",
#' "UseTransactions", "FinalDemand", "InternationalTradeAdjustment,
#' "DomesticUseTransactions", "DomesticFinalDemand",
#' "CommodityOutput, "IndustryOutput", and "DomesticUsewithTrade".
#' @return A list of two-region IO data of model iolevel and year.
getTwoRegionIOData <- function(model, dataname) {
  # Define state, year and iolevel
  if(!"US-DC" %in% model$specs$ModelRegionAcronyms) {
    state <- state.name[state.abb == gsub(".*-", "", model$specs$ModelRegionAcronyms[1])]
  } else {
    state <- "District of Columbia"
  }
  # Define data file name
  filename <- paste(lapply(c("TwoRegion", model$specs$BaseIOLevel, dataname,
                              model$specs$DisaggregationSpecs, model$specs$IOYear,
                              model$specs$IODataVersion),  
                  function(x) x[!is.na(x)]), collapse = "_")
  # Adjust filename to fit what is on the Data Commons
  if(dataname %in% c("UseTransactions", "FinalDemand")) {
    filename <- gsub(dataname, "Use", filename)
  } else if(dataname %in% c("DomesticUseTransactions", "DomesticFinalDemand")) {
    filename <- gsub(dataname, "DomesticUse", filename)
  } else if(dataname %in% c("DomesticUseTransactionswithTrade")){
    filename <- gsub(dataname, "DomesticUsewithTrade", filename)
  } else if(dataname %in% c("UseTransactionswithTrade")){
    filename <- gsub(dataname, "UsewithTrade", filename)
  }
  # Load data
  TwoRegionIOData <- readRDS(loadDataCommonsfile(paste0("stateio/", filename, ".rds")))
  # Keep SoI and RoUS only
  TwoRegionIOData <- TwoRegionIOData[[state]]
  if(dataname %in% c("UseTransactions", "DomesticUseTransactions")) {
    TwoRegionIOData <- TwoRegionIOData[, !(colnames(TwoRegionIOData) 
                                           %in% model$FinalDemandMeta$Code_Loc)]
  } else if(dataname %in% c("FinalDemand", "DomesticFinalDemand")) {
    TwoRegionIOData <- TwoRegionIOData[, model$FinalDemandMeta$Code_Loc]
  } else if(dataname == "ValueAdded") {
    TwoRegionIOData <- TwoRegionIOData[model$ValueAddedMeta$Code_Loc, ]
  }
  return(TwoRegionIOData)
}

#' @description Disaggregate CPI table to ensure the correct dimensions
#' @param df, CPI table
#' @param model An EEIO form USEEIO model object with model specs and IO meta data loaded.
#' @return An expanded CPI table with values replicated for disaggregated sectors.
disaggregateCPI <- function(df, model){
  ## Find rows in IndustryCPI not in IndustryOutput, and duplicate them 
  sector_index <- !(rownames(df) %in% names(model$IndustryOutput))
  disagg_sectors <- rownames(df)[sector_index]
  
  numNewSectors <- (length(model$IndustryOutput) - nrow(df)) / 2
  for (row in disagg_sectors){
    originalIndex <- which(rownames(df)==row)
    originalRowVector <- df[originalIndex,]
    disaggRows <-originalRowVector[rep(seq_len(nrow(originalRowVector)), numNewSectors + 1),,drop=FALSE]
    
    df <- rbind(df[1:originalIndex-1,,drop=FALSE],  #from 1st row to row right before disaggregation
                disaggRows,
                df[-(1:originalIndex),,drop=FALSE])
  }
  return(df)
}


#' @description Generate direct requirements Use table for 2 region models using domestic 
#' Use table with trade data generated by stateior
#' @param model An EEIO form USEEIO model object with model specs and IO meta data loaded.
#' @param domestic A logical parameter indicating whether to DR or Domestic DR.
#' @return A 2-region direct requirements table generated using the domestic Use table with trade
generate2RDirectRequirementsfromUseWithTrade <- function(model, domestic){
  # This function contains code adapted from stateior's validateTwoRegionLagainstOutput()
  # function adapted to work within the useeior package.
  iolevel <- model$specs$BaseIOLevel
  ioschema <- model$specs$BaseIOSchema
  year <- model$specs$IOYear
  state_abb <- sub(".*-","",model$specs$ModelRegionAcronyms[1]) ## Extract characters after -
  
  # Define industries and commodities
  industries <- getVectorOfCodes(ioschema, iolevel, "Industry")
  commodities <- getVectorOfCodes(ioschema, iolevel, "Commodity")
  ita_column <- ifelse(iolevel == "Detail", "F05100", "F051")
  
  if(domestic) {
    ls <- model$DomesticUseTransactionswithTrade
    name <- "Domestic Use table"
  } else {
    ls <- model$UseTransactionswithTrade
    name <- "Use table"
  }
  TwoRegionIndustryOutput <- model$IndustryOutput
  SoI_Industry_Output <- TwoRegionIndustryOutput[endsWith(names(TwoRegionIndustryOutput),
                                                          state_abb)]
  RoUS_Industry_Output <- TwoRegionIndustryOutput[endsWith(names(TwoRegionIndustryOutput),
                                                           "RoUS")]
  # If industry/comm output == 0, it's not viable to generate A matrix, hence set it to 1.
  SoI_Industry_Output[SoI_Industry_Output == 0] <- 1
  
  logging::loginfo(paste0("Generating A matrix of SoI2SoI ", name, " ..."))
  SoI2SoI_A <- normalizeIOTransactions(ls[["SoI2SoI"]][, industries],
                                       SoI_Industry_Output)
  
  logging::loginfo(paste0("Generating A matrix of RoUS2SoI ", name, " ..."))
  RoUS2SoI_A <- normalizeIOTransactions(ls[["RoUS2SoI"]][, industries],
                                        SoI_Industry_Output)
  
  logging::loginfo(paste0("Generating A matrix of SoI2RoUS ", name, " ..."))
  SoI2RoUS_A <- normalizeIOTransactions(ls[["SoI2RoUS"]][, industries],
                                        RoUS_Industry_Output)
  
  logging::loginfo(paste0("Generating A matrix of RoUS2RoUS ", name, " ..."))
  RoUS2RoUS_A <- normalizeIOTransactions(ls[["RoUS2RoUS"]][, industries],
                                         RoUS_Industry_Output)

  U_n_w_trade <- cbind(rbind(SoI2SoI_A, RoUS2SoI_A), rbind(SoI2RoUS_A, RoUS2RoUS_A))
  colnames(U_n_w_trade) <- colnames(model$UseTransactions)
  rownames(U_n_w_trade) <- rownames(model$UseTransactions)
  return(U_n_w_trade)
}


#' Prepares a production demand vector representing production for two region models
#' Demand for SoI = SoI2SoI + RoUS2SoI
#' Demand for RoUS = SoI2RoUS + RoUS2RoUS
#' @param model An EEIO model object with model specs and IO tables loaded
#' @param location, str of location code for demand vector
#' @param domestic A logical parameter indicating whether to generate domestic demand vector.
#' @param demand_type A str indicating whether demand is Production or Consumption
#' @return A named vector with demand
prepare2RDemand <- function(model, location, domestic, demand_type = "Production") {
  # Get state abbreviations, e.g., "US-ME" and "RoUS"
  state_abb <- sub(".*/","",model$FinalDemandMeta$Code_Loc) ## Extract characters after /
  state_abb <- unique(state_abb)
  iolevel <- model$specs$BaseIOLevel
  
  if(domestic){
    use_table <- model$DomesticUseTransactionswithTrade    
  } else {
    use_table <- model$UseTransactionswithTrade
  }
  # Getting list of final demand columns used for the appropriate demand
  if(demand_type == "Production") {
    FD_columns <- unlist(sapply(list("HouseholdDemand", "InvestmentDemand", 
                                     "ChangeInventories", "Export", "Import",
                                     "GovernmentDemand"),
                                getVectorOfCodes, ioschema = model$specs$BaseIOSchema,
                                iolevel = iolevel))
    FD_columns <- FD_columns[FD_columns %in% gsub("/.*", "", colnames(model$FinalDemand))]
    # Calculate production demand for both regions
    ita_column <- ifelse(iolevel == "Detail", "F05100", "F051")
    if(location == state_abb[1]) {
      # calculate production final demand for SoI
      if(domestic) {
        SoI2SoI_y <- rowSums(use_table[["SoI2SoI"]][, c(FD_columns, ita_column, "ExportResidual")])
      } else {
        SoI2SoI_y <- rowSums(use_table[["SoI2SoI"]][, c(FD_columns, "ExportResidual")])
      }
      RoUS2SoI_y  <- rowSums(use_table[["RoUS2SoI"]][, c(FD_columns, ita_column)])
      y_p <- c(SoI2SoI_y, RoUS2SoI_y)
      
    } else if(location == state_abb[2]) {
      # calculate production final demand for RoUS
      if(domestic) {
        RoUS2RoUS_y <- rowSums(use_table[["RoUS2RoUS"]][, c(FD_columns, ita_column, "ExportResidual")])
      } else {
        RoUS2RoUS_y <- rowSums(use_table[["RoUS2RoUS"]][, c(FD_columns, "ExportResidual")])
      }
      SoI2RoUS_y <- rowSums(use_table[["SoI2RoUS"]][, c(FD_columns, ita_column)])
      y_p <- c(SoI2RoUS_y, RoUS2RoUS_y)
    }
    
  } else if(demand_type == "Consumption") {
    # Includes only household, investment, and government consumption as per Ingwersen et al. 2022 (USEEIOv2.0 paper)
    FD_columns <- unlist(sapply(list("HouseholdDemand", "InvestmentDemand", "GovernmentDemand"),
                                getVectorOfCodes, ioschema = model$specs$BaseIOSchema, iolevel = iolevel))
    
    # Calculate consumption demand for both regions
    if(location == state_abb[1]) {
      # calculate production final demand for SoI
      SoI2SoI_y   <- rowSums(use_table[["SoI2SoI"]][, c(FD_columns, "ExportResidual")])
      RoUS2SoI_y  <- rowSums(use_table[["RoUS2SoI"]][, c(FD_columns)])
      y_p <- c(SoI2SoI_y, RoUS2SoI_y)
    } else if(location == state_abb[2]) {
      # calculate production final demand for RoUS
      SoI2RoUS_y  <- rowSums(use_table[["SoI2RoUS"]][, c(FD_columns)])
      RoUS2RoUS_y <- rowSums(use_table[["RoUS2RoUS"]][, c(FD_columns, "ExportResidual")])
      y_p <- c(SoI2RoUS_y, RoUS2RoUS_y)
    }
  }
  names(y_p) <- model$Commodities$Code_Loc
  return(y_p)
}


#' Run validation checks for 2R models and print to console
#' @param model A complete 2R EEIO model: a list with USEEIO model components and attributes
#' @return A list with 2R model results. 
#' @export
print2RValidationResults <- function(model) {
  
  # Check that Production demand can be run without errors
  cat("Checking that production demand vectors do not produce errors for 2-R models, ",
      "as well as validating model components.\n\n")
  printValidationResults(model)
  cat("\n")
  
  # Creating 2-R Production Complete demand vector
  f <- model$DemandVectors$vectors[endsWith(names(model$DemandVectors$vectors), "Production_Complete")][[1]]
  f <- (f + model$DemandVectors$vectors[endsWith(names(model$DemandVectors$vectors), "Production_Complete")][[2]])
  
  cat("Calculating direct results using Production Complete final demand...\n")
  directResultsProductionComplete <- calculateEEIOModel(model, perspective = "DIRECT",
                                                        demand = f, location = NULL,
                                                        use_domestic_requirements = FALSE)
  
  cat("\nCalculating final results using Production Complete final demand...\n\n")
  finalResultsProductionComplete <- calculateEEIOModel(model, perspective = "FINAL",
                                                       demand = f, location = NULL,
                                                       use_domestic_requirements = FALSE)
  
  # Check that Consumption demand can be run without errors
  cat("\n\nChecking that consumption demand vectors do not produce errors for 2-R models.\n\n")
  
  # Creating 2-R Consumption Complete demand vector
  f <- model$DemandVectors$vectors[endsWith(names(model$DemandVectors$vectors), "Consumption_Complete")][[1]]
  f <- (f + model$DemandVectors$vectors[endsWith(names(model$DemandVectors$vectors), "Consumption_Complete")][[2]])
  
  cat("Calculating direct results using Consumption Complete final demand...\n")
  directResultsConsumptionComplete <- calculateEEIOModel(model, perspective = "DIRECT",
                                                         demand = f, location = NULL,
                                                         use_domestic_requirements = FALSE)
  
  cat("\nCalculating final results using Consumption Complete final demand...\n")
  finalResultsConsumptionComplete <- calculateEEIOModel(model, perspective = "FINAL",
                                                        demand = f, location = NULL,
                                                        use_domestic_requirements = FALSE)
  
  twoRegionResults_ls <- list()
  twoRegionResults_ls$directResultsProductionComplete <- directResultsProductionComplete
  twoRegionResults_ls$directResultsConsumptionComplete <- directResultsConsumptionComplete
  twoRegionResults_ls$finalResultsProductionComplete <- finalResultsProductionComplete
  twoRegionResults_ls$finalResultsConsumptionComplete <- finalResultsConsumptionComplete

  # return(twoRegionResults_ls)
}


#' Validate commodity totals between 2R Use table, Make table, and total commodity output objects
#' @param model A complete 2R EEIO model: a list with USEEIO model components and attributes
#' @return A list containing failures of commodity total comparisons between various model objects. 
#' @export
validate2RCommodityTotals <- function(model) {
  
  failures_ls <- list()

  cat("Comparing commodity totals summed from Make and Use (with trade) tables.\n")
  commodityNum <- dim(model$Commodities)[1] # Get number of commodities
  q_make <- colSums(model$V)
  q_use <- rowSums(model$U[1:commodityNum,])#excluding VA rows, including all columns

  failures_ls$Make_Use <- compare2RCommodityTotals(q_make, q_use)
  
  cat("Comparing commodity totals summed from Make and Domestic Use (with trade) tables.\n")
  q_d_use <- rowSums(model$U_d[1:commodityNum,])#excluding VA rows, including all columns
  failures_ls$Make_DUse <- compare2RCommodityTotals(q_make, q_d_use)
  
  
  cat("Comparing commodity totals summed from Make and commodityTotal (model$q) object imported from stateior.\n\n")
  failures_ls$Make_modelq <- compare2RCommodityTotals(q_make, model$q)

  return(failures_ls)
  
}

#' Compare commodity totals between the specified 2R model objects
#' @param q_One A vector of commodity totals derived from specific model object
#' @param q_Two A vector of commodity totals dervied from a different model object than q_One
#' @return A list of sectors that failed the comparison between the two specified q vectors. 
#' @export
compare2RCommodityTotals <- function(q_One, q_Two) {
  
  # Calculate relative differences in q_One and q_Two
  rel_diff_q <- (q_Two - q_One)/q_One
  # Validate relative diff
  validationResults <- formatValidationResult(rel_diff_q, abs_diff = TRUE, tolerance = 0.01)
  failures <- validationResults$Failure
  failuresIndex <- which(rownames(validationResults$RelativeDifference) %in% failures$rownames)
  failures <- cbind(failures, validationResults$RelativeDifference[failuresIndex,1])
  colnames(failures)[1:3] <- c("Commodity", "Validation", "Relative Diff")
  cat(paste(c("Number of failures: ",length(failures$Commodity)),"\n", collapse = " "))
  cat(paste(c("Failing commodities: ", failures$Commodity),"\n", collapse = " "))
  cat("\n")
  return(failures)
  
}