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
  alias <- ifelse(!is.na(model$specs$Alias), model$specs$Alias, NULL)
  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,
                              alias, 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
  logging::loginfo(paste0("Loading ", filename))
  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)
}

#' 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)
}


#' 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 <- unique(model$Industries$Code)
  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) {
    # TODO: CHANGE domestic FROM BOOLEAN TO STRING WITH VALUES 'domestic', 'production', 
    # and 'import', so we can calculate the import matrix in the following if else if else block
    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("/.*", "", model$FinalDemandMeta$Code_Loc)]
    # 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 consumption final demand for SoI
      SoI2SoI_y   <- rowSums(use_table[["SoI2SoI"]][, c(FD_columns)])
      RoUS2SoI_y  <- rowSums(use_table[["RoUS2SoI"]][, c(FD_columns)])
      y_p <- c(SoI2SoI_y, RoUS2SoI_y)
    } else if(location == state_abb[2]) {
      # calculate consumption final demand for RoUS
      SoI2RoUS_y  <- rowSums(use_table[["SoI2RoUS"]][, c(FD_columns)])
      RoUS2RoUS_y <- rowSums(use_table[["RoUS2RoUS"]][, c(FD_columns)])
      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
print2RValidationResults <- function(model) {
  
  # Check that Production demand can be run without errors
  cat("\nChecking that production demand vectors do not produce errors for 2-R models.\n")
  
  if(is.null(model$B)) {
    # Stop validation as no satellite tables
    return()
  }
  
  # 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. 
validate2RCommodityTotals <- function(model) {
  
  failures_ls <- list()

  cat("\nComparing 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 <- compare2RVectorTotals(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 <- compare2RVectorTotals(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 <- compare2RVectorTotals(q_make, model$q)

  return(failures_ls)
  
}

#' Compare  totals between the specified 2R model vectors
#' @param v_One A vector of totals derived from specific 2R model object
#' @param v_Two A vector of totals dervied from a different 2R model object than v_One
#' @return A list of sectors that failed the comparison between the two specified q vectors. 
compare2RVectorTotals <- function(v_One, v_Two) {
  
  # Calculate relative differences in v_One and v_Two
  rel_diff_q <- (v_Two - v_One)/v_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)
  
}


#' Create Make, Use, and Env ratio files for each state from Proxy data for the relevant sectors.
#' @param model An stateior model object with model specs and specific IO tables loaded
#' @param disagg Specifications for disaggregating the current Table
#' @param disaggYear Integer specifying the state model year
#' @param disaggState A string value that indicates the state model being disaggregated. For national models, string should be "US"
#' @return model
createDisaggFilesFromProxyData <- function(model, disagg, disaggYear, disaggState){
  
  # Note: this function assumes: 
  # 1) The disaggregation will use the same proxy values for all disaggregated sectors across all rows and columns. 
  #    That is, if we are disaggregating Summary 22 into the 3 Detail utility sectors, and the proxy allocations are (for example) 0.5/0.25/0.25, then 
  #    in the Use table, the three Detail utility commodities (rows) will have that same split for across all columns (industries/final demand)
  # 2) The disagg parameter will contain a disagg$stateDF variable that includes the data for the relevant disaggState and disaggYear parameters.
  
  if(!is.null(spec$stateFile)){ 
    
    stop("This section of code is meant to be used with 2R models with disaggregated utilities
                           and is not yet fully implemented.")
  }
  
  #Get subset of ratios for current year
  stateDFYear <- subset(disagg$stateDF, Year == disaggYear & State == disaggState)
  
  # If the state/year combination is not found, assume a uniform split between sectors
  if(dim(stateDFYear)[1] == 0){
    
    activity <- unlist(disagg$NewSectorCodes)
    uniformAllocationVector <- 1/length(disagg$NewSectorCodes)
    share <- rep(uniformAllocationVector,length(disagg$NewSectorCodes))
    
    stateDFYear <- data.frame(State = rep(disaggState, length(disagg$NewSectorCodes)),
                              Activity = activity,
                              Share = share,
                              Year = rep(disaggYear, length(disagg$NewSectorCodes)))
    
  }
  
  print(paste0("For ",disaggState,"-",disaggYear, " the allocation to disaggregate ", 
               disagg$OriginalSectorCode, " into ", disagg$NewSectorCodes, " is ", stateDFYear$Share))
  
  # Default Make DF based on proxy employment values
  # Specifying commodity disaggregation (column splits) for Make DF
  industries <- c(rep(disagg$OriginalSectorCode,length(disagg$NewSectorCodes)))
  commodities <- unlist(disagg$NewSectorCodes)
  PercentMake <- stateDFYear$Share # need to add code to ensure that the order of stateDF$Share is the same as the order of disagg$NewSectorCodes
  note <- c(rep("CommodityDisagg", length(disagg$NewSectorCodes)))
  
  makeDF <- data.frame(cbind(data.frame(industries), data.frame(commodities), data.frame(PercentMake), data.frame(note))) #need to rename the columns with the correct column names
  colnames(makeDF) <- c("IndustryCode","CommodityCode",	"PercentMake",	"Note")
  
  
  # Default Use DF based on employment ratios
  # Specifying industry disaggregation (column splits) for Use DF
  industries <- unlist(disagg$NewSectorCodes)
  commodities <- c(rep(disagg$OriginalSectorCode,length(disagg$NewSectorCodes)))
  PercentUse <- stateDFYear$Share
  note <- c(rep("IndustryDisagg", length(disagg$NewSectorCodes)))
  
  useDF <- data.frame(cbind(data.frame(industries), data.frame(commodities), data.frame(PercentUse), data.frame(note))) #need to rename the columns with the correct column names
  useDF_2 <- makeDF # so that colnames match
  colnames(useDF) <- c("IndustryCode","CommodityCode",	"PercentUse",	"Note")
  colnames(useDF_2) <- c("IndustryCode","CommodityCode",	"PercentUse",	"Note")
  
  useDF <- rbind(useDF, useDF_2) #need to bid makeDF because disaggregation procedure requires the UseDF to have the default commodity and industry output.
  
  # Add new DFs to disagg and to model
  disagg$MakeFileDF <- makeDF
  disagg$UseFileDF <- useDF
  
  model$DisaggregationSpecs[[disagg$OriginalSectorCode]] <- disagg
  
  return(model)
  
}